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CHAPTER  1 
INTRODUCTION 

1.1     Motivations 

Noise,  artifacts,  and  low  contrast  can  cause  signal  and  image  degradations  during 
data  acquisition  of  many  signals  and  images,  especially  in  areas  of  medical  image  ap- 
plications. Different  image  modalities  exhibit  distinct  types  of  degradation.  Mammo- 
graphic  images  often  exhibit  low  contrast  while  images  formed  with  coherent  energy, 
such  as  ultrasound,  suffer  from  speckle  noise.  Transmitted  audio  signals  sometimes 
have  the  problem  of  channel-added  noise.  These  degradations  not  only  lower  audio 
or  visual  quality,  but  also  cause  analysis  and  recognition  algorithms  to  sometimes  fail 
to  achieve  their  objectives. 

Since  poor  image  quality  often  makes  feature  extraction,  analysis,  recognition, 
and  quantitative  measurements  problematic  and  unreliable  in  various  areas  of  sig- 
nal/image processing  and  computer  vision,  much  research  has  been  devoted  to  im- 
prove the  quality  of  acquired  signals  and  images  degraded  by  these  factors  [30].  Data 
restoration  techniques  are  often  targeted  to  reduce  noise.  Unlike  noise,  artifacts 
sometimes  comprise  not  only  high  frequency  noise  components,  but  also  middle-to- 
low  frequency  components  which  are  very  hard  to  differentiate  from  other  typical 
features  in  the  spectrum  of  frequency  contents.  Simple  de-noising  techniques  will 
have  problems  reducing  artifacts,  so  application-specific  methods  often  have  to  be 
used  to  eliminate  artifacts  based  on  certain  prior  knowledge.  The  quality  of  low 
contrast  images  can  be  improved  through  designated  feature  enhancement. 


The  promise  of  wavelet  representations  under  dyadic  or  further  sub-octave  wavelet 
transforms  and  the  need  for  improving  degraded  signals/images  have  motivated  this 
dissertation  research.  Reducing  artifacts  is  not  a  major  concern  of  this  research.  The 
major  task  of  this  research  is  to  develop  reliable  techniques  for  reducing  noise  and 
enhancing  salient  features  important  to  various  application  problems.  The  promising 
de-noising  and  feature  enhancement  techniques  may  improve  the  reliability  and  per- 
formance of  signal/image  processing  and  computer  vision  algorithms  for  high-level 
tasks,  such  as  object  detection  or  visual  perception.  In  addition,  we  also  show  the 
capability  of  wavelet  representations  for  other  medical  imaging  applications. 

1.2     Review  of  Related  Methods 

Signal/image  restoration  and  enhancement  have  been  the  focus  of  much  research 
in  the  areas  of  signal  and  image  processing  as  well  as  computer  vision.  Several 
de-noising  methods  and  image  enhancement  techniques  have  been  developed  and  re- 
ported in  the  literature  [30,  15,  26,  17,  62,  58,  68].  Most  of  these  methods  can  be 
roughly  classified  as  spatial-,  statistical-,  and  Fourier-domain-based.  Many  tradi- 
tional methods  were  reviewed  by  Jain  [30].  These  conventional  techniques  for  image 
restoration  and  enhancement  have  shown  certain  limitations  of  balancing  the  effect 
of  removing  noise  and  enhancing  features. 

For  de-noising  purposes,  spatial  and  frequency  domain  smoothing  methods  often 
not  only  reduce  noise,  but  also  smooth  out  high  frequency  components  of  wide-band 
features  as  a  side-effect.  This  is  because  the  smoothing  effect  applies  both  noise  and 
high  frequency  components  of  features.  Over-smoothing  noise  often  causes  certain 
interested  features  being  blurred.  Recently  multiscale  and/or  multiresolution  wavelet 
techniques  for  signal/image  restoration  have  been  reported,  suggesting  improved  re- 
sults in  signal/image  quality  [50,  46,  19,  9,  18,  59,  8,  65]. 


Mallat  and  Hwang  [50]  introduced  a  local  maxima  based  method  for  removing 
white  noise.  Their  method  analyzes  the  evolution  of  local  maxima  of  the  dyadic 
wavelet  transform  cross  scales  and  identifies  the  local  maxima  curves  above  a  cer- 
tain measurement  metric  which  more  likely  correspond  to  features  than  noise.  The 
de-noised  signal/image  is  then  reconstructed  based  on  the  extracted  local  maxima 
corresponding  to  significant  feature  singularity  points.  Lu  et  al.  [46]  further  extended 
the  ideal  of  local  maxima  curves  to  the  local  maxima  tree  cross  scales  and  employed 
a  different  measurement  metric  to  detect  features  from  noise.  Coifman  and  Majid 
[9]  developed  a  wavelet  packet-based  method  for  de-noising  signals.  It  is  an  iterative 
method  for  extracting  features  based  on  the  best  wavelet  packet  basis,  which  removes 
noise  energy  below  a  certain  threshold. 

Donoho  and  Johnstone  [19,  20]  presented  thresholding-based  wavelet  shrinkage 
methods  for  noise  reduction.  These  methods  uniformly  reduce  noise  coefficients  below 
a  global  threshold.  Hard  thresholding  and  soft  thresholding  have  trade  off  between 
preserving  features  and  achieving  smoothness.  Soft  thresholding  de-noising  was  fur- 
ther explained  by  Donoho  [18]  and  proved  that  with  high  probability  the  de-noised 
signal  is  at  least  as  smooth  as  the  noise-free  original  in  a  wide  variety  of  smoothness 
measures  and  comes  nearly  as  close  (in  the  mean  square  sense)  to  the  original  as 
any  other  estimated  results.  But  the  method  still  faces  the  problem  of  balancing  the 
removal  of  noise  and  signal  details  in  order  to  get  a  better  performance  in  terms  of 
visual  quality  and  quantitative  measurements.  Thresholding-based  wavelet  shrinkage 
under  an  orthonormal  wavelet  transform  has  shown  undesired  side-effects,  including 
pseudo-Gibbs  phenomena  [8].  In  order  to  solve  the  problem,  Coifman  and  Donoho 
[8]  developed  a  translation-invariant  de-noising  method.  Their  method  alleviates  the 
problem,  but  oscillations  after  de-noising  remain  visible.   Wei  and  Burrus  [65]  used 


redundant  wavelet  representations  to  achieve  translation-invariant  effects  for  image 
restoration  when  applied  to  various  image  coding  schemes. 

Most  noise  in  reality  is  additive,  but  certain  noise  can  not  be  characterized  well 
as  additive  noise,  such  as  speckle  noise.  Speckle  noise  can  be  better  approximated 
as  multiplicative  noise  [30].  Image  formation  under  coherent  waves  results  in  a  gran- 
ular pattern  known  as  speckle.  The  granular  pattern  is  correlated  with  the  surface 
roughness  of  an  object  being  imaged.  Goodman  [25]  presented  an  analysis  of  speckle 
properties  under  coherent  irradiance,  such  as  laser  and  ultrasound.  The  primary 
differences  between  laser  and  ultrasound  speckle  were  pointed  out  by  Abbott  and 
Thurstone  [1]  in  terms  of  coherent  interference  and  speckle  production.  For  speckle 
reduction,  earlier  techniques  include  temporal  averaging  [25,  1],  median  filtering, 
and  homomorphic  Wiener  filtering  [30].  Homomorphic  Wiener  filtering  is  a  method 
which  converts  multiplicative  noise  into  additive  noise  and  applies  Wiener  low-pass 
filtering  to  reduce  noise.  Similar  to  temporal  averaging,  one  speckle  reduction  tech- 
nique [57]  used  frequency  and/or  angle  diversity  to  generate  multiple  uncorrelated 
synthetic-aperture  radar  (SAR)  images  which  were  summed  incoherently  to  reduce 
speckle.  Hokland  and  Taxt  [28]  reported  a  technique  which  decomposed  a  coherent 
image  into  three  components,  one  of  which,  called  subresolvable  quasi-periodic  scat- 
ter, causes  speckle  noise.  The  component  was  eliminated  by  harmonic  analysis  and 
processing. 

In  the  last  few  years,  several  wavelet  based  techniques  were  developed  for  speckle 
reduction.  Moulin  [54]  introduced  an  algorithm  based  on  the  maximum-likelihood 
principle  and  a  wavelet  regularization  procedure  on  the  logarithm  of  a  radar  image 
to  reduce  speckle.  Guo  et  al.  [27]  first  reported  a  method  based  on  wavelet  shrinkage 
for  speckle  reduction.  In  the  method  of  Guo  et  al,  wavelet  shrinkage  of  a  logarith- 
mically transformed  image  is  applied  for  speckle  reduction  of  SAR  images.    They 


also  proposed  several  approaches  to  combine  data  from  polarization  to  achieve  better 
performance.  We  [71,  72]  have  developed  a  method  for  speckle  reduction  similar  to 
the  one  by  Guo  et  al.  [27].  The  differences  are  that  (a)  noise  is  modeled  as  multi- 
plicative, taking  a  homomorphic  approach  to  reduce  the  noise,  (b)  different  wavelets 
and  multiscale  overcomplete  representations  are  employed  in  our  approach,  and  (c) 
an  enhancement  mechanism  is  incorporated  into  our  de-noising  process.  Thus,  our 
method  can  not  only  reduce  speckle  noise,  but  also  enhance  interesting  features. 

In  the  last  two  decades,  many  image  enhancement  methods  have  been  published 
in  the  literature.  Several  spatial  and  frequency-based  techniques  [11,  30,  24,  44,  58] 
were  developed  for  various  image  enhancement  purposes.  Contrast  stretching,  high- 
pass  filtering,  histogram  modification  methods  are  described  in  Jain  [30].  Contrast 
stretching  was  an  earlier  technique  for  contrast  enhancement  [30].  This  method  has 
limitations  of  selecting  features  based  on  local  information  for  enhancement  because 
it  is  a  global  approach  and  the  enhancement  function  is  linear  or  piece-wise  linear. 
Contrast  stretching  may  also  amplify  noise  when  input  data  are  corrupted  by  noise. 
Some  image  enhancement  schemes  applied  to  medical  image  modalities  have  been 
developed  and  studied  in  the  literature  [30,  58,  45,  39,  35].  Specifically,  spatial  and 
frequency-based  techniques  [30,  24,  11,  44]  have  been  developed  for  ultrasound  image 
enhancement.  A  statistical  enhancement  method,  which  used  the  local  standard 
deviation  of  a  surrounding  region  centered  around  each  pixel  to  replace  its  value  to 
enhance  edges  in  ultrasound  images,  was  reported  by  Geiser  [23]. 

Conventional  filtering-based  techniques  for  de-noising  and  image  enhancement 
have  shown  certain  limited  ability  for  removing  noise  without  blurring  features  and 
for  enhancing  contrast  without  amplifying  noise  because  spatial  and  frequency  rep- 
resentations can  not  separate  features  from  noise  well.  In  the  last  few  years,  many 
techniques  based  on  multiscale  features,  such  as  edges  and  object  boundaries,  have 


achieved  success  for  image  enhancement  in  several  application  areas  [32,  40,  41,  45]. 
Recently  wavelet-based  nonlinear  enhancement  techniques  have  produced  superior 
results  in  medical  image  enhancement  [39,  35,  73]. 

1.3     Objectives  of  De-Noising  and  Enhancement 

To  improve  the  quality  of  acquired  signals  and  images  degraded  by  noise  and/or 
low  contrast,  most  traditional  methods  try  either  to  reduce  noise  or  to  enhance  fea- 
tures. At  first  glance,  de-noising  and  feature  enhancement  appear  to  be  two  conflict- 
ing objectives,  especially  to  traditional  methods  for  de-noising  or  image  enhancement. 
However,  they  are  simply  two  sides  of  the  same  coin.  The  purpose  of  de-noising  is  to 
eliminate  noise,  primarily  in  high  frequency,  while  methods  of  feature  enhancement 
attempt  to  enhance  specific  signal  details,  including  contrast  enhancement.  The  dif- 
ference lies  in  the  fact  that  features  often  occupy  a  wider  frequency  band  than  noise. 
It  is  even  more  difficult  to  achieve  both  objectives  when  feature  details  are  corrupted 
by  noise.  Traditional  spatial  and  filtering-based  methods  for  de-noising  often  reduce 
noise  at  a  price  of  blurring  features  while  single-scale  conventional  methods  for  con- 
trast enhancement  may  amplify  noise.  The  single-scale  representation  of  a  signal  in 
time  (or  pure  frequency)  is  problematic  when  attempting  to  discriminate  signal  from 
noise. 

Because  of  the  limited  ability  of  traditional  techniques  for  de-noising  or  feature 
enhancement,  the  two  conflicting  objectives  can  not  be  accomplished  simultaneously 
through  earlier  methods  under  spatial  or  Fourier  representations  with  a  single  res- 
olution of  frequency  contents.  When  the  two  mechanisms,  de-noising  and  feature 
enhancement,  are  combined  under  a  framework  of  a  new  representation  or  platform 
which  helps  to  overcome  the  drawback  of  each  mechanism  when  acting  alone,  we  will 
have  a  much  better  chance  to  fulfill  the  two  objectives  at  the  same  time.   Wavelet 


transforms  and  wavelet  theory  can  be  one  method  for  new  representation  and  plat- 
form. This  may  be  why  wavelet  representations  have  attracted  more  and  more  at- 
tention of  researchers  aiming  at  signal/image  restoration  and  feature  enhancement. 
Recently  wavelet  based  methods  have  shown  promise  in  accomplishing  the  two  ob- 
jectives at  the  same  time  because  wavelet  decomposition  can  fine  tune  frequency 
resolution  of  signal  details.  We  are  able  to  treat  distinct  components  of  details  at 
fine-to-coarse  scales  differently  to  achieve  desired  effects  of  de-noising  and  feature  en- 
hancement. Algorithms  have  been  developed  under  such  a  multiscale  wavelet  analysis 
framework  [71,  73]. 

1.4     Wavelet  Based  Approaches  for  De-Noising  and  Enhancement 

Since  Morlet  and  Grossmann  [52]  formulated  the  first  wavelet  decomposition, 
wavelet  theory  [52,  12,  13,  14,  6,  47,  48,  49,  51,  64]  has  been  developed  and  well 
documented  in  the  last  14  years.  Some  practical  applications  of  the  theory  have  been 
developed,  but  more  applications  are  still  under  the  developing  stage.  There  are  many 
choices  of  wavelets  with  different  properties  [12].  De-noising  using  some  wavelets  hav- 
ing oscillations  may  lead  to  certain  unwanted  and  undesired  side-effects,  for  example 
noise-induced  ripples  and  oscillations  when  reconstructed  under  incomplete  coeffi- 
cients in  the  wavelet  domain.  This  could  be  one  of  the  major  factors  resulting  in 
artifacts,  including  the  so-called  pseudo-Gibbs  phenomena  in  the  neighborhood  of 
sharp  variation  points  (singularities)  after  de-noising  (for  details  see  Coifman  and 
Donoho  [8]). 

Orthonormal  Wavelet  Transform  (OWT)  and  discrete  Dyadic  Wavelet  Transform 
(DWT)  have  been  used  in  various  application  areas,  such  as  data  compression,  edge 
detection,  texture  analysis,  noise  reduction,  and  image1  enhancement.   The  compact 


and  local  support  of  wavelets  in  spatial  and  frequency  domains  has  been  a  valu- 
able property  for  characterizing  features  locally.  This  enable  us  to  remove  noise  and 
enhance  features  locally  without  affecting  other  features  distant  apart.  In  a  prelim- 
inary implemented  method,  DWT  has  been  adopted  as  our  major  analysis  tool  for 
de-noising  and  contrast  enhancement  [73].  The  reasons  are  quite  obvious.  A  DWT 
with  the  first  order  derivative  of  a  smoothing  function  as  its  basis  wavelet  can  sep- 
arate noise  energy  from  feature  energy  reasonably  well  in  the  wavelet  domain.  The 
DWT  also  correlates  prominent  features  in  its  multiscale  representation,  such  as  edges 
and  object  boundaries.  After  experimental  analysis  and  understanding  of  signal  and 
noise  behaviors  in  scale  space,  we  are  able  to  find  out  which  wavelet  coefficients  to 
modify  to  enhance  certain  features  of  interest  (FOI)  based  on  simple  thresholding 
and  nonlinear  processing.  The  mother  wavelet  is  a  smoothing  function  and  is  anti- 
symmetric with  few  oscillations,  which  keeps  us  relatively  free  from  the  side-effects 
shown  under  OWT  with  a  basis  wavelet  having  slight  oscillations  itself.  This  effect 
can  be  clearly  seen  from  the  de-noised  results  under  OWT  and  our  de-noised  results 
[73].  The  filters  used  to  perform  the  DWT  have  compact  support  of  a  few  taps.  The 
DWT  is  a  stable  and  overcomplete  representation.  DWT  wavelet  coefficients  (WC) 
have  a  more  clear  meaning  that  they  are  proportional  to  the  signal  magnitude  or 
image  intensity  changes  (gradients)  at  certain  scales.  WCs  reflect  energy  in  a  signal, 
so  we  can  rephrase  that  a  DWT  with  the  aforementioned  wavelets  is  a  process  for 
the  diffusion  of  the  energy  of  a  signal  and  converting  it  into  the  energy  of  the  signal 
at  different  scales  in  its  wavelet  representation. 

Even  through  the  DWT-based  algorithm  [73]  has  produced  better  results  than 
de-noising  only  methods  for  signal/image  restoration,  we  have  observed  that  a  DWT 
has  limited  ability  to  characterize  band-limited  features,  such  as  texture  information, 


speech  or  sound  signals  including  ultrasound  signals.  To  more  reliably  identify  fea- 
tures through  the  time-scale  space,  we  formulate  and  implement  a  sub-octave  wavelet 
transform,  which  is  a  generalization  of  the  DWT.  The  sub-octave  wavelet  transform 
provides  a  means  to  visualize  signal  details  in  equal-divided  sub-octave  frequency 
bands  and  is  shown  to  characterize  signal  details  more  effectively.  The  theoretical 
development  of  a  sub-octave  wavelet  transform,  FIR  filter  design,  and  efficient  im- 
plementation are  part  of  this  thesis  research.  Further  more,  in  this  thesis,  we  are 
developing  a  complete  algorithm  and  quantitatively  measure  its  performance  which 
will  be  compared  to  other  techniques  for  de-noising  and/or  feature  enhancement. 

In  an  approach  developed  during  this  research,  we  achieve  de-noising  and  fea- 
ture enhancement  under  a  framework  of  multiscale  sub-octave  wavelet  analysis  and 
judicious  nonlinear  processing  [42,  43].  We  seek  to  eliminate  noise  while  restoring 
or  enhancing  salient  features.  Through  multiscale  representation  by  a  discrete  sub- 
octave  wavelet  transform  (SWT)  with  first  and  second  order  derivative  approxima- 
tions of  a  smoothing  function  as  its  basis  wavelets,  we  can  distinguish  feature  energy 
from  noise  energy  reasonably  well.  The  objectives  of  de-noising  and  feature  enhance- 
ment are  achieved  by  simultaneously  lowering  noise  energy  and  raising  feature  energy 
through  designated  nonlinear  processing  of  wavelet  coefficients  in  the  transform  do- 
main. Through  parameterized  processing,  we  are  able  to  achieve  a  flexible  control  and 
the  potential  to  reduce  speckle  and  restore  (or  even  enhance)  contrast  along  features, 
such  as  object  boundaries.  As  shown  in  our  earlier  work  [73],  this  approach  for  speckle 
reduction  and  contrast  enhancement  is  less  affected  by  pseudo-Gibbs  phenomena  [8]. 

1.5     Organization  of  This  Dissertation 

The  rest  of  this  dissertation  is  organized  as  follows.  In  Chapter  2,  we  review  some 
de-noising  and  enhancement  methods.  We  describe  how  to  regulate  threshold-based 
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wavelet  shrinkage  through  scale  space  and  show  how  to  design  an  enhancement  func- 
tion with  noise  suppression.  In  Chapter  3,  we  present  dyadic  wavelet  transform  based 
techniques  for  de-noising  and  feature  enhancement.  Sample  application  results  and 
analysis  are  presented.  In  Chapter  4,  we  derive  and  formulate  a  sub-octave  wavelet 
transform  mathematically  and  show  how  it  generalizes  the  dyadic  wavelet  transform. 
The  advantage  of  a  sub-octave  representation  over  a  dyadic  wavelet  representation 
is  presented.  Sample  application  results  are  presented.  In  Chapter  5,  we  describe 
how  to  quantitatively  measure  the  performance  of  an  algorithm  for  de-noising  and 
enhancement.  Some  comparisons  are  made  between  the  results  of  other  published 
methods,  and  our  DWT-based  as  well  as  SWT-based  methods.  In  Chapter  6,  we 
apply  wavelet  representations  under  dyadic  and  sub-octave  wavelet  transforms  to 
other  problems  of  medical  image  processing.  Experimental  results  and  analysis  are 
presented.  This  dissertation  is  concluded  in  Chapter  7.  In  Appendix  A,  we  present 
FIR  filters  used  for  dyadic  and  sub-octave  wavelet  transforms.  In  Appendix  B,  we 
introduce  procedures  for  a  fast  implementation  of  a  sub-octave  wavelet  transform  in 
one  and  two  dimensions. 


CHAPTER  2 
DE-NOISING  AND  ENHANCEMENT  TECHNIQUES 

In  this  chapter,  we  are  going  to  overview  the  impact  of  noisy  and  low  contrast  sig- 
nals/images, and  review  some  related  methods  for  de-noising  and  enhancement.  We 
also  describe  how  to  regulate  the  threshold-based  de-noising  techniques  and  design  an 
enhancement  function  with  noise  suppression.  We  then  introduce  the  image  restora- 
tion and  enhancement  techniques  employed  in  our  algorithms.  The  reason  that  we 
put  this  chapter  ahead  of  wavelet  representations  described  in  Chapters  3  and  4  is 
that  both  our  DWT  and  SWT  based  algorithms  share  image  restoration  and  enhance- 
ment techniques  introduced  in  this  chapter.  The  advantage  of  this  organization  is 
that  we  avoid  describing  the  de-noising  and  enhancement  operators  repeatedly  when 
presenting  our  DWT  and  SWT  based  algorithms  for  noise  reduction  and  contrast 
enhancement,  so  we  can  simply  refer  to  the  operators  in  this  chapter. 

2.1     Introduction 

Signal  and  image  degradations  by  noise  and  low  contrast  are  frequent  phenomena 
of  signal/image  data  acquisition,  especially  in  medical  imaging.  Image  degradations 
have  a  significant  impact  on  the  performance  of  human  experts  and  computer-assisted 
methods  for  medical  diagnosis.  For  example,  a  human  medical  expert  may  fail  to 
capture  some  important  information  from  a  noisy  and  low  contrast  image  when  per- 
forming medical  diagnosis,  especially  when  exhausted.  A  cardiologist  may  have  to 
perform  border  interpolation  in  order  to  identify  myocardial  borders  when  border 
information  is  incomplete  and  corrupted  by  speckle  noise  and  make  decision  based 
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on  unreliable  information.  Noise  and  low  contrast  make  it  problematic  for  human 
experts  and  computer  algorithms  to  identify  features  of  diagnostic  importance  in 
medical  imaging.  In  addition,  noise  and  low  contrast  often  make  feature  extraction, 
analysis,  and  recognition  algorithms  unreliable,  so  improving  the  quality  of  acquired 
medical  images  becomes  necessary.  Signal/image  restoration  and/or  enhancement 
are  usually  taken  as  the  first  step  of  a  high  level  task  of  image  processing  and  com- 
puter vision.  For  instance,  in  most  image  segmentation  algorithms,  image  smoothing 
is  usually  carried  out  as  the  first  step  (or  preprocessing)  for  segmentation  in  order  to 
reduce  noise  interference  on  the  performance  of  these  algorithms. 

Most  traditional  approaches  for  de-noising  have  a  single-minded  objective  which 
is  to  reduce  noise  while  minimizing  the  smoothing  of  features.  Traditionally,  noise 
is  frequently  not  a  concern  in  feature  enhancement  algorithms.  In  the  combined  ap- 
proach developed  during  this  thesis  research,  we  will  focus  on  two  goals;  (1)  to  remove 
noise  and  (2)  to  enhance  salient  features  to  a  desired  degree.  As  part  of  this  research, 
we  have  implemented  algorithms  for  removing  additive  and  multiplicative  noise  re- 
spectively while  enhancing  prominent  features  at  the  same  time.  These  algorithms 
are  primarily  based  on  wavelet  representations,  wavelet  shrinkage,  and  feature  em- 
phasis. Wavelet  shrinkage  is  a  technique  which  uniformly  reduces  wavelet  coefficients 
through  a  certain  operator,  such  as  hard  thresholding  or  soft  thresholding.  During 
the  process,  small  coefficients,  mainly  attributed  to  noise,  are  usually  removed.  For 
feature  enhancement,  we  revitalize  low  contrast  FOI  through  feature  emphasis  (in- 
creasing the  energy  level  for  each  of  these  features).  When  the  noise  level  in  a  signal 
or  image  is  high,  these  algorithms  are  capable  of  not  only  removing  noise,  but  also 
restoring  features  to  near  their  original  quality  and  even  enhancing  certain  FOI  se- 
lectively. When  the  noise  level  is  low,  such  as  in  a  low  contrast  medical  image,  our 
algorithms  can  enhance  features  with  noise  suppression  to  avoid  amplifying  noise. 
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The  main  ideas  behind  selected  wavelet  shrinkage  and  salient  feature  emphasis 
encapsulate  the  two  fundamental  objectives  of  de-noising  and  feature  enhancement: 

(a)  Remove  noise,  but  not  features,  and 

(b)  Enhance  the  features  of  interest,  but  not  noise. 

These  are  two  conflicting  objectives  in  the  sense  that  both  sharp  features  and  noise  lie 
in  high  frequency  of  the  spectrum.  Noise  is  often  smoothed  out  at  the  price  of  blurred 
features  left  in  a  traditional  de-noising  algorithm.  On  the  other  hand,  enhancing  cer- 
tain FOI  corrupted  by  noise  is  more  likely  to  amplify  noise  in  an  enhancement-only 
technique  without  a  noise  suppression  mechanism  incorporated.  This  prevents  tradi- 
tional algorithms  from  attempting  to  achieve  both  of  the  objectives  simultaneously 
because  noise  and  features  can  not  be  distinguished  well  in  spatial  or  Fourier  repre- 
sentations. In  Fourier  domain,  de-noising  is  usually  carried  out  through  some  type 
of  low-pass  filtering  in  nature,  including  template-based  spatial  averaging.  On  the 
other  hand,  feature  enhancement  is  accomplished  under  a  certain  type  of  high-pass 
filtering.  They  are  in  conflict  with  each  other  when  performed  on  a  single  set  of  data 
represented  either  in  time  or  in  frequency.  From  this  analysis,  it  sounds  like  that  some 
type  of  band-pass  filtering  may  be  a  choice.  But  in  fact  single  band-pass  filtering  at 
a  frequency  band  has  a  very  limited  capability  of  removing  noise  and  enhancing  fea- 
tures. To  achieve  both  the  objectives,  we  need  a  suitable  representation  or  platform 
which  can  separate  features  from  noise  well  and  locally.  Multiscale  wavelet  represen- 
tation developed  by  Mallat  and  Zhong  [51]  has  shown  promise  in  separating  features 
and  noise  through  scale  space.  As  outlined  by  Daubechies  [13]  as  an  example,  we 
formulated  and  implemented  multiscale  sub-octave  wavelet  representation  [42,  43],  a 
generalization  of  the  dyadic  wavelet  representation,  to  further  improve  the  capability 
of  characterizing  features  from  noise.  A  dyadic  wavelet  transform  is  briefly  explained 
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in  Chapter  3  while  a  sub-octave  wavelet  transform  is  formulated  and  described  in 
Chapter  4. 

2.2     Noise  Modeling 

Noise  modeling  is  an  important  part  of  a  noise  reduction  method  and  it  affects 
which  kind  of  techniques  should  be  used  to  reduce  noise.  Efficient  noise  models  can 
make  de-noising  more  effective.  When  a  noise  behavior  is  not  fully  understood  or 
still  can  not  be  completely  explained,  its  accurate  noise  model  is  very  difficult  to 
obtain.  However,  approximate  noise  models,  such  as  speckle  noise  modeling,  may  be 
used  in  such  a  case.  Continuous  noise  modeling  is  of  theoretical  importance  while 
discrete  noise  models  are  more  related  to  practical  signal/image  processing  for  noise 
reduction.  Through  the  sampling  theory,  a  discrete  noise  model  can  be  obtained 
from  sampling  its  corresponding  continuous  noise  model  with  a  sample  rate  (at  least 
Nyquist  rate)  to  avoid  aliasing  effect. 

2.2.1     Additive  Noise  Model 

For  some  signal/image  processing  applications  considered,  such  as  simulated  sig- 
nals and  mammograms,  noise  is  better  approximated  as  an  additive  phenomenon.  In 
general,  additive  noise  can  be  represented  by  the  formula 

/(x)  =  g{x)  +  ifc(x),  (2.1) 

where  g(x)  is  a  desired  unknown  function.  The  function  /(x)  is  a  noisy  observation 
of  g(x),  ?7a (x)  is  additive  noise,  and  x  is  a  vector  of  spatial  locations  or  temporal 
samples.  By  using  vector  notation,  we  unified  the  noise  model  for  1-D,  2-D,  ...,  iV-D 
cases.  For  our  signal/image  processing,  1-D  and  2-D  noise  models  are  what  we  are 
interested  in.    Noise  77a (x)  is  usually  approximated  as  Gaussian  white  noise,  so  it 
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has  zero  mean  (//  =  0)  and  a  noise  level  a,  the  standard  deviation  of  the  Gaussian 
function.  For  1-D  signal  processing,  we  discretize  Equation  (2.1)  as 

f(n)  =  g(n)  +  %(n),  (2.2) 

where  n  G  Z,  f(n)  =  f(nT  +  s)  {g(n)  and  r)a{n)  are  similar),  T  is  a  sampling  period, 
and  s  is  a  sampling  shift.  For  2-D  image  processing,  Equation  (2.1)  is  discretized  as 

f{m,n)  =  g(m,n)  +  7]a{m,n),  (2.3) 

where  (m,n)  £  Z  ,  /(m,n)  =  f(mTx  +  sx,nTy  +  sy)  (g(m,n)  and  77a(m,  n)  are 
similar),  Tx  and  Ty  are  sampling  periods  along  horizontal  and  vertical  directions,  sx 
and  sy  are  sampling  shifts. 

For  an  additive  noise  model,  there  exist  techniques  based  on  mean  squared  error 
or  I1  norm  optimization  to  remove  noise.  Such  techniques  include  Donoho  and  John- 
stone's wavelet  shrinkage  techniques  [19,  20,  18],  Chen  and  Donoho's  basis  pursuit 
de-noising  [5],  Mallat  and  Hwang's  local-maxima-based  method  for  removing  white 
noise  [50],  and  wavelet  packet-based  de-noising  [9,  10]. 

By  incorporating  de-noising  and  feature  enhancement  mechanisms  within  a  frame- 
work of  wavelet  representations  [73,  42],  we  seek  to  reduce  noise  and  enhance  contrast 
without  amplifying  noise.  We  shall  demonstrate  that  subtle  features  barely  seen  or 
invisible  in  a  mammogram,  such  as  microcalcification  clusters,  spicular  lesions,  and 
circular  (arterial)  calcifications,  can  be  enhanced  via  wavelet  representations  and 
judicious  selection  and  modification  of  transform  coefficients.  Since  our  algorithm 
treats  noise  and  features  independently,  superior  results  were  obtained  for  similar 
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data  when  compared  to  existing  algorithms  designed  for  de-noising  or  enhancement 
alone. 

In  our  investigation,  we  studied  hard  thresholding  and  Donoho  and  Johnstone's 
soft  thresholding  wavelet  shrinkage  [19,  18]  for  noise  reduction.  An  advantage  of  soft 
thresholding  is  that  it  can  achieve  smoothness  while  hard  thresholding  better  pre- 
serves features.  In  our  approach  for  accomplishing  de-noising  and  feature  enhance- 
ment, we  take  advantage  of  both  thresholding  methods.  Donoho  and  Johnstone's 
soft  thresholding  method  [19,  18]  was  developed  on  an  orthonormal  wavelet  trans- 
form [12]  primarily  applied  with  Daubechies's  Symmlet  8  basis  wavelet.  These  results 
showed  some  undesired  side-effects,  from  pseudo-Gibbs  phenomena  [8].  By  using  an 
overcomplete  wavelet  representation  and  basis  wavelets  with  fewer  oscillations,  a  re- 
sult relatively  free  from  such  side-effects  after  de-noising  was  observed  experimentally 
on  similar  data  sets.  In  our  algorithm,  we  first  adapt  regularized  soft  thresholding 
wavelet  shrinkage  to  remove  noise  energy  within  the  finer  levels  of  scale  (such  as  levels 
1  and  2).  We  then  apply  to  wavelet  coefficients  within  the  selected  levels  (such  as 
levels  3  and  4)  of  analysis  a  nonlinear  gain  with  hard  thresholding  incorporated  to 
preserve  features  while  removing  small  noise  perturbations. 

2.2.2     Approximate  Speckle  Noise  Model 

Coherent  interfering  cause  speckle  noise.  An  accurate  and  reliable  model  of  the 
noise  is  desirable  for  efficient  speckle  reduction.  But  it  remains  a  difficult  problem. 
An  approximate  speckle  noise  model  is  formulated  below.  Here,  our  primary  objective 
is  to  accomplish  speckle  reduction  for  2-D  digital  echocardiograms,  so  we  formulate 
the  noise  model  directly  in  two  dimensions.  The  formulation  of  a  one-dimensional 
noise  model  is  similar. 
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Since  speckle  noise  is  not  simply  additive,  Jain  [30]  presented  a  general  model  for 
speckle  noise  as 

f(x,  y)  =  9(x,  y)  rjrhix,  y)  +  r)a(x,  y),  (2.4) 

where  g(x,  y)  is  an  unknown  2-D  function,  such  as  a  noise-free  original  image,  to  be 
recovered,  f(x,y)  is  a  noisy  observation  of  g(x,  y),  77T?l(x,?/)  and  rja{x,y)  are  multi- 
plicative and  additive  noise  respectively,  x  and  y  are  the  variables,  such  as  spatial 
locations,  and  (x,y)  G  R2.  Since  the  effect  of  additive  noise  (such  as  sensor  noise) 
with  level  o-a  is  considerably  smaller  than  multiplicative  noise  (coherent  interfering) 
(||//a(:r,  y)||2  <C  \\g(x,y)7]rh(x1y)\\2)  in  echocardiograms,  Equation  (2.4)  can  be  ap- 
proximated by 

f(x,y)  =  g{x,y)r)rh(x,y).  (2.5) 

To  separate  the  noise  from  the  original  image,  we  take  a  logarithmic  transform  on 
both  sides  of  Equation  (2.5), 

log{f{x,  y))  =  log{g{x,  y))  +  log{r}fh{x,  J/)).  (2.6) 

Equation  (2.6)  can  also  be  rewritten  as 

fl{x,y)=gl{x,y)  +  r}i(x,y).  (2.7) 

Now  we  can  approximate  rjl&(x,y)  as  additive  white  noise  and  may  apply  various 
wavelet-based  approaches  for  additive  noise  reduction.  With  uniform  sampling,  we 
obtain  the  discrete  version  of  equation  (2.7)  as 

fl(m,n)  =gl(m,n)  +  r)l&(m,n),  (2.8) 
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where  (ra,  n)  G  Z2,  /'(ra,  n)  =  fl(mTx  +  sx,  nTy  +  sy),  Tz  and  Ty  are  sampling  peri- 
ods along  horizontal  and  vertical  directions,  sx  and  sy  are  sampling  shifts.  Wavelet 
representation  and  wavelet  transforms  will  be  presented  in  the  next  two  chapters.  To 
describe  how  the  de-noising  method  works,  here,  we  only  need  the  fact  that  a  wavelet 
transform  is  a  linear  transformation,  and  we  borrow  its  notation  for  a  wavelet  rep- 
resentation whose  details  are  given  in  the  following  chapters.  The  symbol  W  is 
represented  as  a  wavelet  transform,  W*  as  a  wavelet  coefficient  at  scale  23  (or  level 
j)  and  direction  d  (1  for  horizontal  and  2  for  vertical),  Sj  is  a  scaling  approximation 
at  a  final  level  J .  By  the  properties  of  a  linear  transform,  we  have 

W[/V  n)]  =  W[gl(m,  n)}  +  W[^(m,  n)]  (2.9) 

after  applying  wavelet  transform  on  the  both  sizes  of  Equation  (2.8)  where 

W[fl(m,n)}  =  {(W^[/V,n)])d=1,2il<,<y,  Sj[fl(m,n)}}, 
W[gl{m,n)]  =  {(Wf\Sf(m1n)])d=1A1<j<Ji  Sj[gl{m,n)}}, 
WhiKn)]  =  {(W/IiA^nJDfciAi^^,  SM(m,n)}}. 

Since  noise  lies  in  high  frequency,  it  will  reduce  to  near  zero  after  a  finite  number 
of  repeated  smoothings,  so  Sj[rjlh(m,n)]  — »  0.  In  fact,  at  most  5  level  wavelet  de- 
composition is  needed  to  smooth  out  noise  for  most  noise  reduction  applications  we 
conducted.  This  is  why  we  only  carry  out  speckle  noise  reduction  through  eliminat- 
ing noise  energy  \J\Wf[rfh{m^  n)])2.  For  image  restoration  purposes,  it  is  desirable 
to  recover  W[gl(m,ri)],  the  wavelet  transform  of  a  desired  function  gl(m,n),  from 
W[/'(ra,n)]  by  reducing  W[r)l&(m,n)}  in  the  wavelet  domain.  By  taking  the  inverse 
wavelet  transform,  we  may  be  able  to  recover  gl(m,n)  or  a  close  approximation.  For 
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noise  reduction  and  feature  enhancement,  we  want  to  increase  further  the  sharpness 
of  features  of  interest,  such  as  myocardial  boundaries,  through  nonlinear  stretching 
for  feature  energy  gain  on  signal  details  Wf[gl{m,n)]. 

Jain  showed  a  similar  homomorphic  approach  [30,  pp.  313-316]  for  speckle  reduc- 
tion of  images  with  undeformable  objects  through  temporal  averaging  and  homomor- 
phic Wiener  filtering.  The  motion  and  deformable  nature  of  human  hearts  through 
time  prevents  us  from  getting  the  same  status  of  the  left  ventricle  for  multiple  frames. 
Because  we  treat  noise  and  feature  components  differently,  we  are  able  to  produce 
a  result  that  is  superior  to  de-noising  only  algorithms.  We  show  that  our  algorithm 
is  capable  of  not  only  reducing  noise,  but  also  enhancing  features  of  diagnostic  im- 
portance, such  as  myocardial  boundaries  in  2-D  echocardiograms  obtained  from  the 
parasternal  short-axis  view. 

2.3     Uniform  Wavelet  Shrinkage  Methods  for  De-Noising 

Threshold-based  de-noising  is  a  simple  and  efficient  technique  for  noise  reduction 
when  being  applied  within  a  framework  of  wavelet  representations  which  can  separate 
features  from  noise.  Hard  thresholding  has  long  been  used  as  a  useful  tool,  includ- 
ing de-noising.  Soft  thresholding  wavelet  shrinkage  for  de-noising  was  developed  by 
Donoho  and  Johnstone  [19].  Hard  thresholding  and  soft  thresholding  have  trade  off 
between  preserving  features  and  achieving  smoothness.  When  features  in  the  wavelet 
domain  can  be  clearly  distinguished  from  noise,  applying  hard  thresholding  wavelet 
shrinkage  can  achieve  a  better  result  than  soft  thresholding.  When  it  is  not  the  case 
and  smoothness  is  more  desirable,  soft  thresholding  should  be  the  choice. 


20 

2.3.1  Hard  Thresholding 

A  hard  thresholding  operation  can  be  expressed  as 

u(x)  =  TH(v(x),t)  =  »(x)(|t/(x)|  >  t)+,  (2.10) 

where  t  is  a  threshold  value,  x£D  where  D  is  the  domain  of  v(x),  and  u(x)  is  the 
result  of  hard  thresholding  and  has  the  same  sign  as  v(x)  if  nonzero.  The  meaning 
of  (|f(x)|  >  t)+  is  defined  by  the  expression 

1    if  |v(x)|  >  t 
0    otherwise. 

2.3.2  Soft  Thresholding 

Soft  thresholding  [19,  18]  is  a  nonlinear  operator  and  can  be  described  by 

u(x)  =  Ts(v(x),t)  =  signMx))  (|t/(x)|  -  t)+i  (2.11) 

where  threshold  parameter  t  is  proportional  to  the  noise  level  and  xGD,  the  domain 
of  v(x),  and  u(x)  is  the  result  of  soft  thresholding  and  has  the  same  sign  as  v(x)  if 
nonzero.  The  expression  (|v(x)|  —  t)+  is  defined  as 


(\v(x)\-t). 


\v(x)\-t      [f\v{x)\>t, 
0  otherwise. 
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Soft  Thresholding  vs  Hard  Thresholding 
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Figure  2.1.  Thresholding  methods:  soft  thresholding  and  hard  thresholding. 


The  function  sign(u)  is  defined  as 


sign(v)  = 


1        if  v  >  0, 
-1     if  v  <  0, 

0       otherwise. 


Figure  2.1  shows  a  soft  thresholding  operation  compared  with  hard  thresholding. 

2.4     Enhancement  Techniques 

In  this  section,  we  describe  how  to  design  an  enhancement  function  with  noise 
suppression  incorporated.  Several  choices  of  enhancement  functions  are  presented. 
Analysis  and  discussion  of  the  reasons  for  our  design  philosophy  are  also  included. 
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Nonlinear  Enhancement  Function: 
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Figure  2.2.  A  nonlinear  gain  function  for  feature  enhancement  with  noise  suppression. 


2.4.1     Enhancement  by  a  Nonlinear  Gain  Function 

In  the  design  of  an  enhancement  function,  we  try  to  accomplish  the  two  tasks  of  an 
effective  enhancement;  (a)  enhance  features  selectively  and  efficiently  and  (b)  avoid 
amplifying  noise.  An  enhancement  operator  with  noise  suppression  is  desirable  and 
can  be  a  choice  for  achieving  the  two  aforementioned  tasks.  Since  we  can  not  fulfill 
the  tasks  satisfactorily  in  the  original  or  Fourier  representation  of  a  signal  or  image, 
this  leads  us  to  look  at  other  representations  through  some  kinds  of  transformation. 
Through  our  study  and  experiments,  we  observed  that  dyadic  wavelet  representations 
have  shown  a  great  promise  for  separating  features  from  noise.  Therefore,  we  can 
apply  the  kind  of  enhancement  functions,  which  will  be  introduced  momentarily, 
to  enhance  FOI.  We  thereafter  generalize  dyadic  wavelet  representations  to  produce 
sub-octave  wavelet  representations  for  characterizing  band-limited  features  frequently 
seen  in  medical  images  more  efficiently. 
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A  parameterized  nonlinear  gain  function,  which  is  targeted  to  accomplish  the  two 
tasks  of  an  effective  enhancement,  can  be  formulated  as 


ENL(v)  =  < 


0  if|v|<Ti, 

ssign{v)(T2  +  T{{\v\-T2)/T)a)    if  T2  <  \v\  <  T3,  (2-12) 

s  v  otherwise, 


where  v  £  [—  1,  I],  0  <  a  <  1,  T  =  (T3  —  T2),  s  is  a  positive  scaling  factor  which 
is  used  to  adjust  the  overall  energy  of  a  processed  image.  Parameters  T1;  T2,  and 
T3  are  selected  values.  For  each  input  value  v  less  than  T\,  the  small  coefficient  is 
more  likely  resulted  from  noise  where  v  is  a  normalized  coefficient.  For  input  value 
v  greater  than  T3,  the  contrast  of  the  corresponding  feature  of  v  is  already  relatively 
high.  No  special  treatment  for  the  coefficient  is  needed,  so  we  only  do  linear  scaling 
which  is  needed  to  keep  the  enhancement  function  from  becoming  decreasing,  which 
may  cause  artifacts.  The  normalized  coefficients  within  the  range  between  T2  and 
T3  are  what  we  would  like  to  enhance  because  their  contrast  is  relatively  low  and 
our  features  of  interest  have  the  corresponding  coefficients  in  this  range.  Thus,  we 
nonlinearly  stretch  their  energy  gain  to  revitalize  these  features.  The  range  between 
T\  and  T2  is  considered  as  a  risk  area.  Both  noise  and  features  may  have  components 
in  this  range,  so  we  try  to  avoid  amplifying  noise  by  simply  linear  scaling  and  treated 
this  range  similar  to  the  area  of  values  greater  than  T3.  Figure  2.2  shows  a  sample 
enhancement  function.  This  enhancement  operator  is  less  flexible  than  the  operator 
in  Section  2.4.2,  but  it  is  computationally  more  efficient.  The  operator  can  serve  as 
a  choice  if  speed  is  a  concern. 
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Generalized  Adaptive  Gain:  C=1 0,  B=0.35,  T1  =0.1 ,  T2=0.2,  T3=0.9 


Figure  2.3.  A  generalized  adaptive  gain  function. 
2.4.2     Enhancement  by  Generalized  Adaptive  Gain 

In  the  last  few  years,  several  wavelet-based  enhancement  techniques  have  been  de- 
veloped [46,  45,  37,  39].  Adaptive  gain  through  nonlinear  processing  [37,  34]  has  been 
used  to  enhance  features  in  digital  mammograms.  The  adaptive  gain  through  non- 
linear processing  is  further  generalized  [42]  to  incorporate  hard  thresholding  in  order 
to  avoid  amplifying  noise  and  to  remove  small  noise  perturbations.  The  generalized 
adaptive  gain  (GAG)  nonlinear  operator  is  defined  as 


'GAG 


(v) 


0  if  \v\  <  7\, 

q  [T2  +  a  (sigm(c(u  -  b))  -  sigm(-c(w  +  b)))}    if  T2  <\v\<  T3,  (2.13) 

s  v  otherwise, 


where  v  6  [-1,  1],  a  =  (T3  -T2)a,  q  =  ssign(v),  u  =  (\v\  -  T2)/(T3  -  T2),  b  G  (0,  1), 
0  <  T\  <  T2  <  T3  <  1,  c  is  a  gain  factor,  s  is  a  scaling  factor  which  is  used  to  adjust 
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the  overall  energy  of  a  processed  image.  Parameter  a  can  be  computed  by 


a       sigm(c(l  -  b))  -  sigm(-c(l  +  b)) '  (2'H) 


sigm(^)  =  1  +  e_v-  (2-15) 

Parameters  7\,  T2,  and  T3  are  selected  values.  When  T\  =  T2  =  0  and  T3  =  1, 
the  expression  is  equivalent  to  the  adaptive  gain  nonlinear  function  used  previously 
[37,  34].  The  interval  [T2,  T3]  serves  as  a  sliding  window  for  feature  selectivity.  The 
slide  can  be  adjusted  to  emphasize  coefficients  within  a  specific  range  of  variation. 
To  increase  the  overall  energy  of  a  processed  image,  we  assign  s  a  value  greater  than 
one  (s  >  1).  Similarly  we  may  reduce  the  energy  of  a  processed  image  by  letting 
s  <  1.  When  s  =  1,  the  scaling  factor  s  does  not  contribute  to  the  overall  energy 
change,  and  makes  the  overall  operator  equivalent  to  the  operator  reported  by  Laine 
and  Zong  [42].  Thus,  by  selecting  gain  values  and  local  windows  of  energy,  we  may 
achieve  a  more  desirable  enhancement.  Figure  2.3  presents  an  adaptive  gain  function 
with  feature  selectivity. 


CHAPTER  3 
DYADIC  WAVELET  REPRESENTATIONS 

In  this  chapter,  we  describe  multiscale  wavelet  transforms  at  dyadic  scales  adopted 
in  our  preliminary  algorithms  for  de-noising  and  enhancement  as  well  as  edge  detec- 
tion. Wavelet-based  de-noising  and  enhancement  are  presented  in  this  chapter  while 
edge  detection  through  wavelet  maximum  representation  will  be  described  in  Chap- 
ter 6.  Since  we  use  a  dyadic  wavelet  transform  primarily  for  discrete  signal  and 
digital  image  processing,  the  transformation  is  presented  in  the  discrete  domain  from 
an  implementation  point  of  view.  For  our  purposes,  we  are  only  interested  in  the 
overcomplete  (redundant)  representation  under  a  finite-level  discrete  dyadic  wavelet 
transform.  For  more  theoretical  work  and  continuous  dyadic  wavelet  transforms,  Mal- 
lat  and  Zhong  have  presented  in-depth  details  [51,  69].  DWT-based  algorithms  for 
signal/image  processing  are  developed  on  dyadic  wavelet  representations  described 
in  this  chapter,  image  restoration  and  enhancement  operators  presented  in  the  last 
chapter. 

3.1     Discrete  Dyadic  Wavelet  Transform 

The  discrete  dyadic  wavelet  transform  developed  by  Mallat  and  Zhong  [51]  has 
been  previously  applied  to  areas,  including  edge  detection,  texture  analysis,  noise 
reduction,  and  image  enhancement.  Multiscale  representation  under  a  dyadic  wavelet 
transform  provides  a  useful  framework  for  characterizing  features  in  terms  of  sharp 
variation  points.  For  de-noising  and  enhancement  purposes,  compactly  supported 
wavelets  can  be  utilized  to  eliminate  noise  and  sharpen  contrast  within  structures 
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and  along  object  boundaries  without  affecting  distant  features  [39,  35,  73]  because 
wavelet  transforms  localize  feature  representations. 

3.1.1     One-Dimensional  Dyadic  Wavelet  Transform 

First,  we  describe  the  discrete  dyadic  wavelet  transform  in  one  dimension  and 
then  extend  it  to  two  dimensions.  For  discrete  signal  and  digital  image  processing, 
only  a  finite-level  discrete  dyadic  wavelet  transform  is  usually  needed  for  practical 
applications.  A  J-level  discrete  dyadic  wavelet  transform  of  a  1-D  discrete  function 
f(n)  €  12{Z)  can  be  represented  as 

W[/(n)]  =  {(VWWDi^xj,  Sj[f(n)}},  (3.1) 

where  Wj[/(ri)]  is  a  wavelet  coefficient  at  scale  2J  (or  level  j),  location  uGZ,  Sj[f(n)] 
is  a  coarse  scale  approximation  at  the  final  level  J  and  position  n.  Wavelet  coefficient 
Wj[f(ri)}  and  scaling  approximation  Sj[f(n)]  at  level  j  can  be  defined  as 

+oo 

Wj\f(n)]  =  f*ifo(n)  =    £    /M^(n-n'),  (3.2) 

n'=— oo 

+oo 

Sj[f{n)}  =  f  *  win)  =     £    f(n')<Pv(n  -  n')>  (3-3) 


respectively,  where  ip2j(n)  =  ^I^ifi)  and  V?2j (n)  =  ^^(27)  are  analysis  wavelets 
and  scaling  functions  dilated  at  scale  2J .  The  inverse  discrete  dyadic  wavelet  trans- 
form can  be  represented  as  W_1  (W_1[W[/(n)]]).  For  a  perfect  decomposition  and 
reconstruction,  we  have 


f(n)  =  W-llW[f(n)]}  =  Sj[f(n)}  *  fa{n)  +  £  \Vj[f(n)}  *  j2i(n),  (3.4) 
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where  </v(n)  =  ip2J{~ n)-  ^n  order  to  get  a  perfect  reconstruction  of  a  1-D  discrete 
function,  analysis  and  synthesis  wavelets  and  the  scaling  function  should  satisfy 

l^)|2    =    E^(2M7(2M  +  I^(27u;)|2 

3=1 

+00 

=    ^^(2^)7(2^).  (3.5) 

The  finite-level  dyadic  wavelet  decomposition  in  (3.1)  forms  a  complete  representa- 
tion for  a  J-level  dyadic  wavelet  transform.  For  a  particular  class  of  dyadic  wavelets, 
such  as  the  first  order  derivatives  of  spline  smoothing  functions,  the  finite-level  di- 
rect and  inverse  discrete  dyadic  wavelet  transform  of  a  1-D  discrete  function  can  be 
implemented  in  terms  of  three  filters,  H,  G,  and  K.  The  three  filters  should  satisfy 
the  following  condition 

\H(u)\2  +  G{u;)K{lj)  =  1,  (3.6) 

where  H(u)  is  a  low  pass  filter,  G(co)  and  K(u)  are  high  pass  filters. 

The  dyadic  wavelet  decomposition  in  Equation  (3.1)  can  be  formulated  in  terms 
of  the  following  recursive  relations  between  two  consecutive  levels  j  and  j  +  1  in  the 
Fourier  domain  as 

I  WJ+1[f(u)}  =  G(2M^[/M],  (3-7) 

Sj+1[f(u)]  =  H(Vu;)SJlf(u;)],  (3.8) 

where  j  >  0,  and  S0[f(uj)}  :=  f(uj).  The  reconstruction  So[f(co)]  from  a  dyadic 
wavelet  decomposition  can  be  represented  recursively  as  (j  changes  from  J  —  1  to  0) 

Sjlfiu)}  =  Wj+1lf(u>)]K@u)  +  Sj+1[f(u)]H(2^),  (3.9) 
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Figure  3.1.  A  3-level  DWT  decomposition  and  reconstruction  of  a  1-D  function. 

where  H  is  the  complex  conjugate  of  H .  The  DWT  decomposition  and  reconstruction 
based  on  the  above  recursive  relations  are  shown  as  a  block  diagram  in  Figure  3.1. 
The  process  of  wavelet  decomposition  is  referred  to  as  wavelet  analysis  while  the 
wavelet  reconstruction  process  is  sometimes  called  wavelet  synthesis. 

3.1.2     Two-Dimensional  Dyadic  Wavelet  Transform 

In  the  rest  of  this  section,  we  shall  present  discrete  dyadic  wavelet  analysis  and 
synthesis  of  a  2-D  discrete  function  (image).  The  decomposition  will  produce  both 
high-to-middle  frequency  signal  details  (wavelet  coefficients)  and  a  low  frequency 
scaling  approximation  (scaling  coefficients)  of  an  image  at  some  final  level  of  analysis. 

Similarly  a  finite-level  discrete  dyadic  wavelet  transform  is  desirable  for  our  digital 
image  processing.  A  J-level  discrete  dyadic  wavelet  transform  of  a  2-D  discrete 
function  f(nx,ny)  G  12{Z2)  can  be  represented  as 


W[f(nx,ny)}  =  {(Wyi/^n^JfciAi^j,  Sj[f(nx,  />,,)]}.  (3.10) 
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where  W?[f(nx,ny)]  is  a  wavelet  coefficient  at  scale  23  (or  level  j),  position  {nx,ny), 
and  spatial  orientation  d  (1  for  horizontal  and  2  for  vertical),  Sj[f(nx,  ny)\  is  a  coarse 
scale  approximation  at  the  final  level  J  and  position  {nx,ny).  Wavelet  coefficient 
Wj[f(nx,ny)]  and  scaling  approximation  Sj[f(nx,ny)]  at  level  j  can  be  defined  as 

+oo  +oo 

Wf[f{nx,ny)]  =  f*i/>*j(nx,ny)=     £       £    f(m',n')^(m  -m',n-  n'), (3.11) 

m'  =  —  oo  n'  —— oo 

+oo  +0O 

S'i[/(nI,n2/)]  =  f  *ip2j{nx,ny)  =     ^        ]T    /(m',  n')tp2J{m  -  m',n  -  n'),(3.12) 


m'—— oo  n'  =— oo 


respectively,  where  if)$j{nx,ny)  =  ^J^df ,  ^)  and  </?2j  (n*,  %)  =  ^T^df  >  ^)  are 
analysis  wavelets  and  scaling  functions  dilated  at  scale  2J,  and  d  =  1,2  represents 
horizontal  or  vertical  spatial  orientation.  In  our  approach  for  de-noising  and  enhance- 
ment, we  are  interested  in  some  basis  wavelets  which  are  the  first  order  derivatives 
of  continuous  smoothing  functions  d{x, y);  thus,  ipl(nx,ny)  and  ip2(nx,ny)  can  be 
formulated  as 


V  {nx,ny)  = 


dx 


2  dd(x,y) 

,  ip  {nx,ny)  = 


dy 

y  =  ny  a 


x  =  nx 
y  =  ny 


(3.13) 


Convolution  with  dilated  il>l(nx,ny)  and  ip2{nx,ny)  produces  sharp  variations  along 
horizontal  and  vertical  directions  for  salient  features.  In  wavelet  frame  represen- 
tations, we  can  employ  a  different  synthesis  basis  wavelet  7d(nx,ny)  for  the  recon- 
struction of  the  original  2-D  discrete  function.  The  inverse  discrete  dyadic  wavelet 
transform  can  be  represented  as  W_1  (W~l[W[f(nx,ny)}]).  For  a  perfect  decompo- 
sition and  reconstruction, 


f(nx,ny)    =    W^lWlfin^ny)]] 

J       2 

=    Sj[f{nx,ny)\  *ip2->{nx,ny)  +  J2  Y,  wfU(n^  ny)\  *7&(n*>n2/)>  (3-14) 

3=1  d=l 


31 

where  ifi2J{nx,  ny)  —  (P2J{~nxi  ~ny)-  m  order  to  get  a  perfect  reconstruction  of  a  2-D 
discrete  function,  analysis  and  synthesis  wavelets,  the  scaling  function  should  satisfy 

J      2 

\${u>x,uy) |2  =  EE  i^d{^ujx^uy)Y{2]uJx^ujy)  +  |0(2Jw„2Ju;,)|2.         (3.15) 

The  finite-level  dyadic  wavelet  decomposition  in  (3.10)  generates  a  complete  repre- 
sentation for  a  J-level  dyadic  wavelet  transform.  For  a  particular  class  of  2-D  dyadic 
wavelets,  such  as  the  first  order  directional  derivatives  of  spline  smoothing  functions, 
Mallat  and  Zhong  [51]  showed  that  the  finite-level  direct  and  inverse  dyadic  wavelet 
transform  of  a  2-D  discrete  function  can  be  implemented  in  terms  of  four  1-D  filters, 
H,  G,  K,  and  L.  The  four  filters  should  satisfy  the  following  perfect  decomposition 
and  reconstruction  conditions 

\H(uj)\2  +  G(u)K(uj)  =  1,  (3.16) 

1  +  I//HI2 

L{u)  = .  (3.10 

Mallat  and  Zhong  [51]  also  showed  how  to  design  1-D  finite  impulse  response  (FIR) 
filters,  H,  G,  K,  and  L,  for  a  2-D  wavelet  transform. 

Similar  to  the  1-D  case,  a  2-D  dyadic  wavelet  decomposition  in  Equation  (3.10)  can 
be  formulated  in  terms  of  the  following  recursive  relations  between  two  consecutive 
levels  j  and  j •  +  1  in  the  Fourier  domain  as 

W}+l[f(ux,uy)]  =  GiVuJSAfiu^ujy)],  (3.18) 

W?+1[f(ux,uy)]  =  G(2*«r)Si[/(w,,t*)],  (3.19) 

Sj+1[f(ux,uy)]  =  H{2^x)H(2'ujy)Sj[f(uxyUy)]i  (3.20) 
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Figure  3.2.  A  3-level  DWT  decomposition  and  reconstruction  of  a  2-D  function. 

where  j  >  0,  and  So[f(u;x,ujy)]  :=  f{cux,ujy),  the  Fourier  transform  of  f{nx,ny).  The 
reconstruction  So[f(ux,u;y)]  from  a  dyadic  wavelet  decomposition  can  be  represented 
recursively  as 

Sj[f{ux,uy)]  =  W}+1[f(ux,uy)]K(2jujx)L{2juy)  +  W]+l\f{ujxlujy)}L{2^x)K{2^y) 

+Sj+1lf{ux,uy)]H{Vux)H(yuy),  (3.21) 


where  H  is  again  the  complex  conjugate  of  H.  A  2-D  DWT  decomposition  and 
reconstruction  based  on  the  above  recursive  relations  are  shown  as  a  functional  block 
diagram  in  Figure  3.2  for  J  =  3.  For  a  pair  of  2-D  analysis  and  synthesis  filter  banks 
shown  in  Figures  3.3  and  3.4,  reconstructed  f*(nx,ny)  is  equal  to  f(nx,ny)  when  no 
processing  is  performed  on  W[/(nx,ny)].  The  2-D  analysis  and  synthesis  filter  banks 
in  Figures  3.3  and  3.4  are  constructed  using  FIR  filters  shown  in  Table  3.1  where,  for 
instance,  H{u)  =  £n  h{n)e~lu}n. 
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Figure  3.3.  A  2-D  analysis  filter  bank. 


34 


Profiles 


d=l 


d=2 


J=l 


j=2 


j=3 


1  -.  . 

ft. 

V    _  1 

dc-cap  (j=3) 


Figure  3.4.  A  2-D  synthesis  filter  bank. 
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Table  3.1.  Impulse  responses  of  filters  H(u),  G{uj),  K(u),  and  L(u>). 


FIR  filters  for  m  =  4  and  c  =  2 

n 

h(n) 

g(n) 

k{n) 

l(n) 

-4 

0.001953125 

-3 

-0.00390625 

0.015625 

-2 

0.0625 

-0.03515625 

0.0546875 

-1 

0.25 

1.0 

-0.14453125 

0.109375 

0 

0.375 

-1.0 

-0.36328125 

0.63671875 

1 

0.25 

0.36328125 

0.109375 

2 

0.0625 

0.14453125 

0.0546875 

3 

0.03515625 

0.015625 

4 

0.00390625 

0.001953125 

3.2     DWT-Based  De-Noising  and  Feature  Enhancement 

In  this  section,  we  first  introduce  algorithms  for  noise  reduction  and  feature 
restoration  or  enhancement  based  on  an  additive  noise  model  presented  in  Section 
2.2.1  and  for  speckle  reduction  with  feature  enhancement  based  on  an  approximate 
speckle  noise  model  formulated  in  Section  2.2.2.  We  then  describe  the  methods  and 
present  formulation  for  de-noising  and  enhancement  based  on  the  operators  intro- 
duced in  Sections  2.3  and  2.4. 

3.2.1     Algorithm  for  Additive  Noise  Reduction  and  Enhancement 

Because  de-noising  and  enhancement  techniques  are  incorporated  into  a  frame- 
work of  wavelet  representations  under  dyadic  wavelet  transforms,  our  algorithm  for 
noise  reduction  and  contrast  enhancement  consists  of  four  major  steps.  In  these 
steps,  parameterized  de-noising  and  enhancement  operators  are  utilized.  The  sample 
experimental  results  are  shown  for  these  operations.  The  parameters  can  be  fine 
tuned  to  achieve  two  distinct  purposes.  One  is  for  de-noising  with  feature  restoration 
while  the  other  is  for  image  enhancement  with  noise  suppression.  They  are  related 
in  certain  sense,  such  as  removing  noise  and  improving  the  quality  of  features.  These 
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methods  are  designed  to  remove  noise  with  feature  restoration  or  enhancement  in  an 
additive  noise  model.  The  four  steps  of  a  DWT-based  de-noising  and  enhancement 
method  are  listed  as  follows: 

1.  Carry  out  a  DWT  to  obtain  a  complete  representation  of  noisy  data  in  the 
wavelet  domain. 

2.  Shrink  transform  coefficients  within  the  finer  scales  to  partially  remove  noise. 

3.  Emphasize  features  through  a  nonlinear  point-wise  operator  to  increase  energy 
among  features  within  a  specific  range  of  variation. 

4.  Perform  an  inverse  DWT  and  reconstruct  the  signal/image. 

Unlike  Donoho  and  Johnstone's  methods  [19]  for  de-noising,  an  advantage  of  this 
method  is  that  it  also  applies  feature  enhancement  to  further  improve  the  performance 
of  signal/image  restoration.  This  algorithm  has  an  ability  to  suppress  noise  (without 
amplifying  noise)  when  applied  for  contrast  enhancement  compared  to  enhancement 
only  methods. 

3.2.2     Algorithm  for  Speckle  Reduction  with  Feature  Enhancement 

Speckle  noise  was  modeled  as  approximate  multiplicative  noise  in  Section  2.2.2. 
Similar  to  the  method  in  Jain  [30],  we  apply  a  homomorphic  approach  to  reducing 
speckle  noise.  The  algorithm  consists  of  six  major  steps.  The  six  steps  of  a  DWT- 
based  de-noising  and  enhancement  method  for  the  speckle  noise  model  are  listed  as 
follows: 

1.  Perform  a  logarithmic  transform  to  convert  an  image  containing  multiplicative 
noise  into  an  image  with  additive  noise. 

2.  Carry  out  a  DWT  and  obtain  a  complete  representation  of  the  log-transformed 
image  in  the  transform  domain. 


37 

3.  Shrink  coefficients  within  the  finer  scales  to  partially  remove  noise  energy. 

4.  Emphasize  features  through  nonlinear  point-wise  processing  to  increase  the 
energy  among  features  within  a  specific  range  of  variation. 

5.  Perform  an  inverse  DWT  and  reconstruct  the  de-noised  and  enhanced  image  so 
that  it  approximates  its  noise-free  original  in  log  scale  with  features  enhanced. 

6.  Finally,  perform  an  exponential  transform  on  the  reconstructed  image  to  con- 
vert it  from  log  scale  to  linear  scale.  The  resulting  image  is  now  de-noised  and 
enhanced. 

This  method  takes  a  similar  homomorphic  transform  to  convert  multiplicative  noise 
into  additive  noise.  Unlike  Jain's  method  [30],  we  incorporate  a  feature  enhancement 
mechanism  into  the  noise  reduction  procedure  to  sharpen  blurred  features  (feature 
restoration  or  enhancement)  after  de-noising. 

Wavelet  representations  under  discrete  dyadic  wavelet  transforms  were  described 
in  Section  3.1  for  both  one  and  two  dimensions.  These  de-noising  and  contrast 
enhancement  schemes  are  based  on  wavelet  shrinkage  and  feature  emphasis  on  top 
of  the  wavelet  representations.  Wavelet  shrinkage  is  a  technique  to  uniformly  reduce 
wavelet  coefficients  in  order  to  remove  noise  coefficients  for  the  purpose  of  de-noising. 
Feature  emphasis,  on  the  other  hand,  is  trying  to  increase  the  magnitudes  of  feature's 
coefficients  to  gain  energy  for  low  contrast  features.  Below,  we  describe  how  to 
perform  DWT-based  de-noising  and  enhancement. 

3.2.3     DWT-Based  De-Noising 

Since  dyadic  wavelet  transforms  with  first  order  derivatives  of  smoothing  functions 
as  basis  wavelets  can  efficiently  identify  features  with  sharp  variation,  we  are  able 
to  achieve  the  objective  of  noise  reduction  through  either  hard  thresholding  or  soft 
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thresholding.  Hard  thresholding  preserves  features  better  while  soft  thresholding  can 
achieve  the  effect  of  smoothness. 

Here  we  describe  threshold-based  de-noising  in  two  dimensions.  One  dimensional 
case  is  similar.  To  achieve  the  purpose  of  de-noising  through  hard  thresholding,  we 
can  modify  D WT  coefficients  for  noise  reduction  by 

W*>*[f(nx,  ny)]  =  Mf  TH(Wf[f{nx,  ny)]/Mf,  tj) ,  (3.22) 

Mf  =  max(\Wf\f(nx,n9)]\),  (3.23) 

where  d  =  1, 2,  j  =  1, ...,  k,  k  <  J,  and  fj  is,  in  general,  a  threshold  related  to  noise 
level  and  scale.  Parameter  td  can  be  directionally  related  if  we  have  orientation  prefer- 
ence. Th  is  the  hard  thresholding  operator  presented  in  Section  2.3.1.  The  threshold 
td  should  be  selected  to  possibly  remove  most  noise  coefficients  while  preserving  fea- 
ture coefficients.  Selection  of  thresholds  in  [71,  73]  was  trial-and-error  based.  The 
selection  can  be  guided  by  examining  the  histogram  and  energy  distribution  of  coef- 
ficients. Wavelet  transforms  generate  a  small  number  of  large  coefficients  carrying  a 
significant  amount  of  energy,  especially  from  fine  to  coarse  scales,  for  sharp  variation 
points  while  producing  a  large  number  of  small  coefficients  mostly  corresponding  to 
noise.  Thresholds  decrease  from  fine  to  coarse  scales  because  noise  energy  is  smoothed 
out  through  repeated  smoothings  (scaling)  by  low  pass  filtering.  This  point  is  made 
clear,  as  shown  in  Figure  3.5  for  Donoho  and  Johnstone's  synthetic  signal  "Blocks". 
The  guideline  is  to  select  decreasing  thresholds  which  can  remove  a  great  number  of 
small  coefficients  carrying  most  noise  energy  and  keep  a  limited  number  of  large  coef- 
ficients for  feature  energy.  Thresholds  through  fine  to  coarse  levels  may  be  regulated 
by  a  decreasing  function  which  will  be  discussed  momentarily. 
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Figure  3.5.  Coefficient  and  energy  distributions  of  signal  "Blocks", 
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When  noise  level  is  high,  hard  thresholding  de-noising  may  not  be  able  to  achieve 
overall  smoothness.  If  this  is  the  case,  we  can  carry  out  soft  thresholding  based  de- 
noising.  For  reducing  noise  and  achieving  smoothness  effect,  we  can  modify  DWT 
coefficients  by 

W}>*[f(nx,ny)}  =  MfTs(W?[f(nx,ny)]/Mf),td3) ,  (3.24) 

Mf  =  max{\W*[f(nx,ny)}\),  (3.25) 

J  nx,ny  J 

where  d  =  1,  2,  j  =  1, ...,  k,  k  <  J,  and  td  is  a  threshold  usually  related  to  noise  level 
and  scale.  Ts  is  the  soft  thresholding  operator  presented  in  Section  2.3.2.  Thresholds 
can  be  selected  similarly  based  on  the  above  discussion.  Wavelet  coefficients  are 
normalized  to  the  range  between  -1  to  1  before  thresholding  operations. 

3.2.4     Regulating  Threshold  Selection  through  Scale  Space 

Donoho  and  Johnstone's  method  of  soft  thresholding  uses  a  single  global  threshold 
[19,  20]  under  orthonormal  wavelet  transforms.  Since  noise  coefficients  under  a  DWT 
have  average  decay  through  fine-to-coarse  scales,  we  can  regulate  both  soft  and  hard 
thresholding  by  applying  coefficient  dependent  thresholds  at  different  scales.  The 
regulated  threshold  td  can  be  computed  through  a  linearly  decreasing  function 


td  =  i 


{Tmax  -  oc(j  -  1))  a*    if  Tmax  -  a(j  -  1)  >  Tmin  , 
Tmin  ad  otherwise, 


(3.26) 


where  od  is  the  standard  deviation,  a  is  a  decreasing  factor  between  two  consecutive 
levels,  Tmax  is  a  maximum  factor  related  to  od  while  Tmin  is  a  minimum  factor, 
1  <  j  <  J,  and  d  G  {1,2}.  When  the  noise  level  in  original  corrupted  data  is 
unknown,  some  methods  use  the  standard  deviation  to  approximate  the  noise  level,  so 
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max 


J  Level 

Figure  3.6.  A  sample  scaling  factor  function. 

the  thresholds  are  related  to  aj.  Figure  3.6  shows  a  sample  scaling  factor  function  for 
the  computation  of  regulated  thresholds.  Our  de-noising  algorithms  are  implemented 
in  a  way  that  a,  Tmax,  and  Tmin  can  be  interactively  tuned  to  achieve  distinct  effect 
of  noise  reduction. 

3.2.5     DWT-Based  Enhancement  with  Noise  Suppression 

Through  either  a  nonlinear  gain  function  or  generalized  adaptive  gain  nonlin- 
ear operator,  we  can  achieve  the  effect  of  contrast  enhancement  for  certain  FOI  by 
processing  DWT  coefficients  as 


W}'*[f(nx,ny)]  =  MfE0p(Wf[f(nx,ny)]/Mf), 


(3.27) 


Mf  =  max(\W?[f(nx,ny)]\), 


nx,ny 


(3.28) 


where  position  (nx,ny)  £  D,  the  domain  of  f(nx,ny),  d  =  1,2,  j  €  {k ./}.  and 

1  <  k  <  J.  The  enhancement  operator  Eop  can  be  ENL  or  Eqag  presented  in 
Section  2.4.  Since  these  operators  are  defined  on  input  values  between  -1  to  1,  we 
normalize  wavelet  coefficients  before  applying  the  operators. 
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(b)  Reconstructed  aflcr  SoftThrcsholding 
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Figure  3.7.  Pseudo-Gibbs  phenomena,  (a)  Orthonormal  wavelet  transform  of  an  orig- 
inal signal  and  its  noisy  signal  with  added  spike  noise,  (b)  Pseudo-Gibbs  phenomena 
after  both  hard  thresholding  and  soft  thresholding  de-noising  under  an  OWT. 
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Figure  3.8.  Multiscale  discrete  wavelet  transform  of  an  original  and  noisy  signals. 
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The  Enhanced  Signal  and  the  Processed  DWT  Coefficients 
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Figure  3.9.   DWT-based  reconstruction  after  (a)  hard  thresholding,  (b)  soft  thresh- 
olding, and  (c)  soft  thresholding  with  enhancement. 
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3.3     Application  Samples  and  Analysis 

In  this  section,  we  present  the  experimental  results  of  our  algorithms  on  some 
sample  signals  and  images.  First,  we  show  that  these  algorithms  are  less  affected 
from  pseudo-Gibbs  phenomena  through  a  simple  and  illustrative  experiment.  Our 
results  are  compared  to  Donoho  and  Johnstone's  results  to  show  this  effect.  We 
then  present  de-noising  and  enhancement  results  for  additive  noise  and  speckle  noise 
models. 

3.3.1     Less  Affection  from  Pseudo-Gibbs  Phenomena 

Coifman  and  Donoho  [8]  showed  that  both  hard  and  soft  thresholding  de-noising 
under  an  orthonormal  wavelet  transform  produced  undesired  side-effects,  including 
pseudo-Gibbs  phenomena.  To  solve  the  problem,  they  [8]  presented  translation- 
invariant  de-noising  methods  to  overcome  the  artifacts  partially  caused  by  the  lack 
of  shift-invariance  of  an  OWT.  Their  methods  alleviated  the  problem  by  making  it 
less  obvious,  but  oscillations  after  de-noising  remained  visible.  Several  experimental 
results  showed  that  our  algorithms  were  less  affected  from  pseudo-Gibbs  phenomena 
[73] .  We  have  used  a  simple  and  intuitive  synthetic  signal  to  demonstrate  this  effect  of 
our  algorithms  when  compared  to  Donoho  et  a/.'s  methods.  Our  experimental  results 
on  Donoho  and  Johnstone's  four  synthetic  signals  also  demonstrate  this  point. 

Figures  3.7,  3.8,  and  3.9  are  used  to  show  that  our  methods  are  relatively  free 
from  pseudo-Gibbs  phenomena.  We  generated  a  synthetic  signal  to  illustrate  what 
may  cause  the  side-effect  and  how  our  methods  can  basically  avoid  it.  Figures  3.7(a) 
shows  an  original  signal,  its  noisy  version  with  added  spike  noise,  and  the  orthonor- 
mal wavelet  coefficients  of  the  original  and  noisy  signals.  Figure  3.7(b)  shows  the 
effect  of  pseudo-Gibbs  phenomena  under  Donoho  et  a/.'s  hard  thresholding  and  soft 
thresholding  methods  through  WaveLab  (a  software  package  from  Donoho's  research 
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group).  Notice  that  a  feature  of  sharp  variation  produces  not  only  large  coefficients 
but  also  small  coefficients  under  an  OWT.  The  small  coefficients  are  removed  under 
both  hard  and  soft  thresholding  methods.  A  typical  orthonormal  wavelet  usually  has 
at  least  certain  oscillations  [12,  64]  in  order  to  satisfy  the  (admissibility)  condition  of 
an  orthonormal  wavelet 

I^MI2 


/  duo  <  oo, 

J-oo  \u)\ 


where  ip(uj)  =  -4=  f*™  ^(x)elxuJdx.  In  the  spatial  domain,  the  corresponding  contin- 
uous wavelet  function  ip{x)  has  sufficient  decay  and  satisfy 


/+oo 
ip{x)dx  =  0. 
-oo 


Figures  3.9(a)  and  3.9(b)  present  our  de-noising  results  under  regulated  hard 
thresholding  and  soft  thresholding.  Both  methods  remove  the  noise  without  causing 
pseudo-Gibbs  phenomena,  but  soft  thresholding  also  smoothes  features  (step  edges) 
a  little  bit.  The  features  are  basically  restored  through  our  enhancement  mechanism 
in  figure  3.9(c).  This  experiment  is  used  to  illustrate  the  fact  that  the  results  of  our 
algorithm  are  less  affected  from  the  side-effect  (pseudo-Gibbs  phenomena)  compared 
to  the  results  from  Donoho  et  a/.'s  methods. 

3.3.2     Additive  Noise  Reduction  and  Enhancement 

Based  on  an  additive  noise  model,  here,  we  present  application  results  of  de- 
noising  and  enhancement  for  synthetic  signals  and  medical  images.  The  first  part  is 
targeted  for  signal/image  restoration  while  the  second  part  is  for  enhancement  with 
noise  suppression. 
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Figure  3.10.  De-Noised  and  feature  restored  results  of  DWT-based  algorithms;  first 
row:  original  signal,  second  row:  noisy  version,  third  row:  de-noised  only  result,  and 
fourth  row:  de-noised  and  enhanced  result  signal. 

De-Noising  with  Feature  Restoration 

In  this  part  of  experiments,  our  de-noising  and  enhancement  techniques  are  pri- 
marily used  for  signal/image  restoration.  To  achieve  this  objective,  we  want  to  reduce 
noise  and  restore  salient  features.  Since  noise  reduction  usually  causes  features  to  be 
blurred,  our  enhancement  methods  are  deployed  to  sharpen  the  blurred  features. 

Figures  3.10  displays  de-noised  with  feature  restored  results  of  our  DWT-based 
algorithms.  First  rows  of  figure  3.10(a)-(d)  are  original  signals,  second  rows  are  noisy 
signals,  third  rows  are  de-noised  only  results,  and  the  last  rows  are  de-noised  with 
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Figure  3.11.  De-Noising  and  enhancement,  (a)  Original  signal,  (b)  Signal  (a)  with 
added  noise  of  2. 52dB.  (c)  Soft  thresholding  de-noising  only  (11.07dB).  (d)  De-Noising 
with  enhancement  (12.25dB). 


(a) 


(b) 


(c) 


(d) 


Figure  3.12.  De-Noising  and  enhancement,  (a)  Original  image,  (b)  Image  (a)  with 
added  noise  of  2.5dB.  (c)  Soft  thresholding  de-noising  only  (11.75dB).  (d)  De-Noising 
with  enhancement  (14.87dB). 
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feature  restored  result  signals.  Our  results  are  pretty  close  to  Coifman  and  Donoho's 
new  and  best  results  produced  by  the  full  cycle-spinning  translation-invariant  de- 
noising  algorithm  which  is  computationally  more  complex.  For  their  "Bumps"  and 
"Doppler"  test  signals,  our  results  are  better  than  their  best  results.  When  processing 
a  noised  signal  or  image  with  low-contrast,  this  algorithm  can  be  used  to  boost  up 
the  contrast  through  adding  external  energy  to  its  signal  energy  by  specifying  a  larger 
gain  factor. 

Figures  3.11  and  3.12  show  our  de-noising  and  enhancement  results  on  Mallat 
and  Hwang's  test  signal  and  image.  The  signal  and  image  are  corrupted  by  noise 
of  higher  levels,  but  we  get  better  results  (higher  recovery  SNR)  than  Mallat  and 
Hwang's  results  on  the  same  signal  and  image  with  less  noise. 

Enhancement  with  Noise  Suppression 

In  this  part  of  experiments,  we  try  to  achieve  contrast  enhancement  without 
amplifying  noise.  Figures  3.13  and  3.14  show  the  de-noising  and  enhancement  results 
on  two  MRI  head  images  with  unknown  noise  level.  The  experimental  results  of  de- 
noised  and  enhanced  images  from  our  algorithm  are  visibly  and  quantitatively  better 
than  the  results  from  the  thresholding-based  methods  alone  for  de-noising,  especially 
for  high  level  noise. 

3.3.3     Speckle  Reduction  with  Feature  Enhancement 

Speckle  reduction  and  contrast  enhancement  can  be  accomplished  in  the  trans- 
form domain  by  judicious  multiscale  nonlinear  processing  of  wavelet  coefficients 
{Wf[f(nx,ny)])d=i^,i<j<j  obtained  under  dyadic  wavelet  transforms.  Through  the 
approximate  speckle  noise  model  in  Section  2.2.2,  we  can  usually  separate  the  noise 
component  from  a  desired  function.  Wavelet  transforms  help  to  further  distinguish 
signal  from  noise  in  the  spatial-scale  space. 


49 


(a) 


(b) 


(c) 


Figure  3.13.  De-Noising  and  enhancement,  (a)  Original  MRI  image,  (b)  De-Noising 
only,  (c)  DWT-based  de-noising  with  enhancement. 


(a) 


(b) 


(c) 


Figure  3.14.  De-Noising  and  enhancement,  (a)  Original  MRI  image,  (b)  De-Noising 
only,  (c)  DWT-based  de-noising  with  enhancement. 


Figure  3.15.  An  algorithm  for  speckle  reduction  and  contrast  enhancement. 
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Figure  3.16.  Results  of  de-noising  and  enhancement,  (a)  A  noisy  ED  frame,  (b) 
Wavelet  shrinkage  de-noising  only  method,  (c)  DWT-based  de-noising  and  enhance- 
ment. 

Our  multiscale  homomorphic  algorithm,  as  shown  in  Figure  3.15,  for  speckle  re- 
duction and  feature  enhancement  was  tested  on  echocardiograms  of  varying  quality. 
These  image  sequences  were  acquired  from  the  parasternal  short-axis  view.  Figures 
3.16  and  3.17  show  the  results  of  de-noising  with  or  without  feature  enhancement 
on  end  diastolic  (ED)  and  end  systolic  (ES)  frames.  The  speckled  original  frames 
are  shown  first.  Results  from  wavelet  shrinkage  de-noising  only  and  de-noising  with 
enhancement  are  shown  in  the  Figures  3.16(b)  and  3.16(c)  respectively.  Figure  3.18 
shows  a  nonlinear  operator  for  enhancing  the  image  in  Figure  3.17(a).  This  operator 
looks  much  different  from  Figure  2.3  because  of  the  log  transform  effect.  Experimen- 
tal results  are  also  compared  with  other  speckle  reduction  techniques,  such  as  median 
filtering,  extreme  sharpening  combined  with  median  filtering  [11,  44],  homomorphic 
Wiener  filtering,  and  a  wavelet  shrinkage  de-noising  [19,  18].  Figures  3.19  and  3.20 
show  sample  results  of  the  above  mentioned  methods  on  two  typical  frames  from  two 
different  echocardiographic  sequences  with  the  left  ventricle  as  the  region  of  interest. 
Figure  3.19(a)  is  sample  noisy  image.  The  result  of  median  filtering  with  a  5x5 
mask  is  shown  in  Figure  3.19(b).  Figure  3.19(c)  displays  sample  result  of  extreme 
sharpening  combined  with  median  filtering.    The  result  from  homomorphic  Wiener 
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(a) 


(b) 


(c) 


Figure  3.17.  Results  of  de-noising  and  enhancement,  (a)  A  noisy  ES  frame,  (b) 
Wavelet  shrinkage  de-noising  only  method,  (c)  DWT-based  de-noising  and  enhance- 
ment. 

filtering  is  shown  in  Figure  3.19(d).  The  last  two  images  in  Figures  3.19(e)  and 
3.19(f),  display  the  results  from  wavelet  shrinkage  de-noising  only  and  our  de-noising 
and  enhancement  algorithms.  The  algorithm  produces  smoothness  inside  a  uniform 
region  and  improves  contrast  along  structure  and  object  boundaries  in  addition  to 
speckle  reduction.  The  de-noised  and  enhanced  results  of  noisy  echocardiographic 
images  from  this  algorithm  appear  to  outperform  the  results  from  soft  thresholding 
de-noising  alone.  Our  current  algorithm  is  implemented  such  that  most  parameters 
are  set  dynamically  for  adaptive  de-noising  and  feature  enhancement. 

3.4     Clinical  Data  Processing  Study 

A  study  of  clinical  image  processing  was  conducted  to  investigate  the  effect  of 
de-noising  on  the  consistency  and  reliability  to  manually  defined  borders  of  the  left 
ventricle  in  2-D  short-axis  echocardiographic  images  [70].  Experimental  results  in- 
dicate the  algorithm  is  promising.  Myocardial  borders  manually  defined  by  expert 
observers  exhibit  less  variation  after  de-noising.  It  seems  that  in  echocardiograms, 
where  no  real  borders  are  clearly  visible  and  incomplete,  expert  borders  usually  in- 
dicate a  close  range  where  real  borders  may  occur.   When  two  expert  borders  agree 
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Figure  3.18.  A  generalized  adaptive  gain  function  for  processing  an  echocardiogram 
in  Figure  3.17(a). 

with  each  other,  the  range  of  real  borders  is  more  likely  limited  around  the  two  expert 
borders.  The  study  of  clinical  image  processing  shows  that  de-noising  and  feature 
enhancement  help  to  improve  the  consistency  and  reliability  of  manually  defined 
borders  by  expert  observers. 

The  set  of  test  images  in  our  study  of  clinical  image  processing  was  selected 
from  an  echocardiographic  database  exhibiting  diverse  image  quality.  Sixty  systolic 
sequences  of  2-D  short-axis  echocardiographic  images  were  selected.  Half  of  the  test 
images  were  rated  as  good  quality  while  the  rest  were  considered  as  poor  quality. 
For  more  details  about  how  these  echocardiographic  sequences  were  acquired,  we 
refer  the  reader  to  Wilson  and  Geiser  [66].  Statistical  results  have  shown  that  there 
is  some  improvement  in  consistency  and  reliability  for  manually  defined  borders  by 
expert  observers  examining  de-noised  images  compared  to  their  original  noisy  images. 
Quantitative  measurements  were  calculated  in  terms  of  the  mean  of  absolute  border 
differences  (MDistDiff)  in  distance  (mm)  and  the  mean  of  border  area  differences 
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(b) 


(c) 


(d) 


(e) 


(f) 


Figure  3.19.  Results  of  various  de-noising  methods,  (a)  Original  image  with  speckle 
noise,  (b)  Median  filtering,  (c)  Extreme  sharpening  combined  with  median  filtering, 
(d)  Homomorphic  Wiener  filtering,  (e)  Wavelet  shrinkage  de-noising  only,  (f)  DWT- 
based  de-noising  with  enhancement. 
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(a) 


(b) 


(c) 


(d) 


(e) 


(f) 


Figure  3.20.  Results  of  various  de-noising  methods,  (a)  Original  image  with  speckle 
noise,  (b)  Median  filtering,  (c)  Extreme  sharpening  combined  with  median  filtering, 
(d)  Homomorphic  Wiener  filtering,  (e)  Wavelet  shrinkage  de-noising  only  method. 
(f)  DWT-based  de-noising  and  enhancement. 
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Table  3.2.  Quantitative  measurements  of  manually  defined  borders. 


All  Test  Images 

Good  Images 

Poor  Images 

Ori     vs     Enh 

Ori     vs     Enh 

Ori     vs     Enh 

MDistDiff    Endo  (in  mm) 

2.1040      1.8168 

1.5972      1.5322 

2.6118      2.1014 

Epi     (in  mm) 

1.7846      1.6743 

1.3979      1.5886 

2.1713      1.7601 

MAreaDiff  Endo  (in  cm2) 

2.3731      1.8893 

1.6597      1.4543 

3.0865      2.2058 

Epi     (in  cm2) 

2.5676      2.0799 

1.5823      1.9540 

3.5528      2.3243 

(MAreaDiff)  in  cm2.  The  border  difference  was  measured  by  its  close  approximation 
in  64  radial  directional  difference  from  an  estimated  center  [66]  of  the  left  ventricle. 
Manually  defined  borders  by  experts  looking  at  poor  images  were  improved  more 
than  those  of  good  images  after  de-noising.  The  statistical  results  of  quantitative 
measurements  of  two  sets  of  expert  manually  defined  borders  are  shown  in  Table  3.2. 
The  statistical  computation  results  listed  under  the  column  "Ori"  are  the  quantitative 
measurements  between  two  sets  of  expert  borders  on  the  original  speckled  images 
while  the  results  under  the  column  "Enh"  are  based  on  the  de-noised  and  enhanced 
images.  It  is  worth  mentioning  that  a  single  set  of  de-noising  and  enhancement 
parameters  were  used  to  process  all  the  test  echocardiographic  images  used  in  this 
study.  We  suggest  that  a  single  value  set  of  parameters  should  be  enough  for  de- 
noising  and  enhancing  a  class  of  images  with  a  similar  noise  pattern  and  selected 
features. 

Figure  3.21  shows  the  correlation  between  the  areas  delineated  by  the  two  expert 
observers.  The  four  diagrams  in  Figure  3.21(a)  present  the  correlation  of  ED  Epi 
(epicardial)  border  areas,  ES  Epi  border  areas,  ED  Endo  (endocardial)  border  areas, 
and  ES  Endo  border  areas  on  the  original  noisy  images.  The  four  diagrams  in  Fig- 
ure 3.21(b)  show  similar  results  for  the  de-noised  images  with  features  restored  or 
enhanced.  The  solid  lines  in  the  figure  are  the  linear  regression  lines,  while  the  dash 
and  dotted  lines  are  ideal  regression  lines.    From  the  diagrams,  it  is  clear  that  the 
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ED  Epi  Area  -  Original  Image 


ES  Epi  Area  -  Original  Image 


Observer  1 
ES  Endo  Area  -  Original  Image 


ED  Epi  Area  -  DcNoiscd  Image 


ES  Epi  Area  -  DcNoiscd  Image 
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Observer  1 

ED  Endo  Area  -  DcNoiscd  Image 


Observer  1 
ES  Endo  Area  -  DcNoiscd  Image 


0  10  20  30 

Observer  I 


(a) 


(b) 


Figure  3.21.  Area  correlation  between  manually  defined  borders  by  two  expert  car- 
diologist observers. 

points  which  represent  the  two  expert  border  areas  on  the  same  de-noised  image  are, 
in  general,  more  toward  the  ideal  regression  line.  Additional  improvement  can  be  seen 
on  the  Endo  area  correlation  for  the  de-noised  images.  For  most  echocardiograms  in 
the  study,  there  is  usually  less  Endo  border  information  than  Epi  border  information. 
Noisy  border  information  (low  signal-to-noise  ratio  (SNR))  affects  border  interpola- 
tion by  human  observers  for  the  manually  defined  borders.  After  de-noising,  Endo 
border  information  in  terms  of  SNR  is  improved,  so  the  expert  border  areas  tend  to 
agree  with  each  other,  especially  ES  Endo  areas.  The  statistical  computation  results 
shown  in  Table  3.2,  show  evidence  for  this  analysis. 

Figure  3.22  shows  the  distributions  of  mean  border  differences  on  the  original 
images;  (a)  the  distribution  of  Epi  ED  border  differences,  (b)  the  distribution  of  Epi 
ES  border  differences,  (c)  the  distribution  of  Endo  ED  border  differences,  and  (d)  the 
distribution  of  Endo  ES  border  differences.  Figures  3.23(a)-(d)  show  the  distributions 
of  mean  border  differences  on  the  enhanced  images,  similar  to  Figures  3.22(a)-(d). 
The  solid  lines  in  Figures  3.22  and  3.23  are  the  third  order  polynomial  fitting  curves  in 
a  least-squares  sense.  With  the  same  scale  for  both  Figures  3.22  and  3.23,  Figure  3.23 
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Epi  ED  Border  Difference  Dlslnbulion  (Graph  Dale   24-Jul-97) 
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Endo  ED  Border  Difference  Distribution  (Graph  Date   24-Jul-97) 
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Epi  ES  Border  Difference  Distribution  (Graph  Dale   24-Jul-97) 
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Endo  ES  Border  Difference  Dlslnbulion  (Graph  Dale:  24-Jul-97) 
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Figure  3.22.  Border  difference  variation  on  the  original  images,  (a)  Distribution 
of  Epi  ED  border  differences,  (b)  Distribution  of  Epi  ES  border  differences,  (c) 
Distribution  of  Endo  ED  border  differences,  (d)  Distribution  of  Endo  ES  border 
differences.  The  solid  lines  are  the  third  order  polynomial  fitting  curves. 
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Table  3.3.    Quantitative  measurements  of  inter-observer  mean  border  differences  in 
mm  on  original  versus  enhanced  images,  as  shown  in  Figure  3.25. 


Epi-ED     Epi-ES 

Endo-ED     Endo-ES 

Ori 

4.7 

4.0 

6.3 

4.6 

Enh 

1.4 

2.3 

1.2 

1.3 

shows  that  border  distance  differences  for  enhanced  images  have  smaller  means  and 
standard  deviations  than  the  corresponding  differences  for  the  original  noisy  images 
as  shown  in  Figure  3.22. 

Figure  3.24  shows  an  example  of  de-noising  and  image  enhancement.  Figure  3.25 
shows  the  same  images  as  Figures  3.24  with  two  expert  manually-defined  borders 
overlaid.  Significant  overall  improvement  on  the  agreement  of  two  expert  borders  is 
visible  from  the  overlaid  borders  of  the  enhanced  images  compared  to  the  original 
images,  as  shown  visually  in  Figure  3.25  and  quantitatively  in  Table  3.3.  The  Endo 
borders  have  more  improvement  than  the  Epi  borders  based  on  quantitative  mea- 
surements and  visual  appearance.  Statistical  analysis  shows  improvement  in  terms  of 
the  mean  of  absolute  border  differences  and  mean  border  area  differences  of  de-noised 
images  compared  to  their  original  images.  From  the  overall  statistical  analysis,  the 
greatest  impact  is  on  the  expert  borders  drawn  on  images  with  poor  image  quality, 
such  as  the  images  in  Figure  3.24. 


59 


Epi  ED  Border  Difference  Distribution  (Graph  Date:  29-Jul-97) 
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Endii  ED  Border  Difference  Dislnbulion  (Graph  Dale:  29-Jul-97> 
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Endo  ES  Border  Difference  Dislnbulion  (Graph  Dale:  29-Jul-97) 
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Figure  3.23.  Border  difference  variation  on  the  enhanced  images,  (a)  Distribution 
of  Epi  ED  border  differences,  (b)  Distribution  of  Epi  ES  border  differences,  (c) 
Distribution  of  Endo  ED  border  differences,  (d)  Distribution  of  Endo  ES  border 
differences.  The  solid  lines  are  the  third  order  polynomial  fitting  curves. 


60 


(a)  Original  ED 


(b)  Original  ES 


(c)  Enhanced  ED 


(d)  Enhanced  ES 


Figure  3.24.  De-noising  and  image  enhancement:   (a)  An  original  ED  frame;  (b)  An 
original  ES  frame;  (c)  The  enhanced  ED  frame;  (c)  The  enhanced  ES  frame. 
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(a)  Original  ED 


(b)  Original  ES 


(c)  Enhanced  ED 


(d)  Enhanced  ES 


Figure  3.25.  Image  and  border  display:  (a)  An  original  ED  frame  with  manually- 
defined  borders  overlaid;  (b)  An  original  ES  frame  with  manually-defined  borders 
overlaid;  (c)  The  enhanced  ED  frame  with  manually-defined  borders  overlaid;  (c) 
The  enhanced  ES  frame  with  manually-defined  borders  overlaid.  Red  and  yellow 
borders  represent  the  two  observers. 


CHAPTER  4 
SUB-OCTAVE  WAVELET  REPRESENTATIONS 

In  this  chapter,  we  introduce  sub-octave  wavelet  transforms.  A  sub-octave  wavelet 
transform  is  a  generalization  of  a  traditional  dyadic  wavelet  transform  and  we  use 
an  example  to  show  the  advantage  of  sub-octave  representations  over  dyadic  wavelet 
representations  for  characterizing  band-limited  features  frequently  seen  in  medical 
imaging.  We  formulate  both  continuous  and  discrete  sub-octave  representations  in 
one  and  two  dimensions. 

4.1     Introduction 

Our  DWT  based  algorithms  for  de-noising  and  enhancement  have  achieved  im- 
proved performance  compared  to  other  published  methods.  Through  our  analysis 
and  experiments,  we  observed  that  a  DWT  has  a  limited  ability  to  characterize  fea- 
tures, such  as  texture,  and  subtle  features  of  importance  in  mammographic  images. 
The  traditional  DWT  is  an  octave-based  transform  where  scales  increase  as  powers 
of  two  [51].  However,  the  best  representation  of  a  signal's  details  may  exist  be- 
tween two  consecutive  levels  of  scale  within  a  DWT  [36].  To  more  reliably  isolate 
noise  and  identify  features  through  scale  space,  we  designed  a  multiscale  sub-octave 
wavelet  transform  (SWT),  which  generalizes  the  DWT.  A  sub-octave  wavelet  trans- 
form provides  a  means  to  represent  details  within  sub-octave  frequency  bands  of 
equally-spaced  divisions  within  each  octave  frequency  band.  The  theoretical  devel- 
opment of  a  sub-octave  wavelet  transform  and  its  efficient  implementation  was  briefly 
described  by  Laine  and  Zong  [42],  and  later  explained  and  extended  [43]. 
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The  initial  motivation  for  us  to  explore  sub-octave  decomposition  was  that  we  had 
observed  the  limitation  of  dyadic  wavelet  transforms  for  characterizing  band-limited 
features  and  sought  better  frequency  resolution  for  detecting  such  subtle  features.  A 
dyadic  wavelet  transform  is  an  octave-based  transformation,  where  scales  increase  as 
powers  of  two.  Daubechies  [13]  introduced  the  generalization  and  extension  of  wavelet 
decomposition  and  reconstruction  under  the  context  of  orthonormal  wavelet  trans- 
forms by  sub-band  splitting  and  presented  early  examples.  An  extension  and  gener- 
alization of  dyadic  wavelet  transforms  is  multiscale  sub-octave  wavelet  decomposition 
and  reconstruction.  Both  Daubechies's  methods  and  our  techniques  for  sub-octave 
wavelet  transforms  have  achieved  similar  (better)  frequency  localization.  However, 
we  are  primarily  interested  in  overcomplete  (redundant)  wavelet  representations,  a 
generalization  of  dyadic  wavelet  representations.  Most  orthonormal  wavelet  bases 
have  the  effect  of  decorrelating  features  while  dyadic  wavelet  bases  correlate  salient 
features  through  scales,  which  is  what  we  are  most  interested  in  for  enhancement 
purposes.  In  the  rest  of  the  chapter,  we  present  the  mathematical  formulation  of 
sub-octave  wavelet  bases  in  both  space  and  frequency  domains.  The  decomposition 
and  reconstruction  procedure  is  carried  out  in  terms  of  filter  bank  theory  and  band 
splitting  techniques. 

4.2     Continuous  Sub-Octave  Wavelet  Transform 

Through  a  wavelet  transform,  a  function  (input  signal)  can  be  represented  by 
its  projection  onto  a  family  of  wavelet  bases  for  decomposition  and  possibly  perfect 
reconstruction.  If  the  family  of  wavelets  bases  {ipn(x)}  is  complete  and  orthonormal, 
a  wavelet  transform  with  critical  sampling  is  usually  referred  to  as  an  orthonormal 
wavelet  transform  [13].  A  Haar  wavelet  is  a  simple  example  of  an  orthonormal  wavelet . 
However,  the  Haar  wavelet  is  a  discontinuous  function  and  is  not  localized  in  the 
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frequency  domain.  The  analysis  filters  {H,  G}  for  computing  an  orthonormal  wavelet 
transform  must  satisfy  the  following  design  constraints 

\H(u)\2  +  \G(u)\2  =  l, 

H(uj)  G{uo)  =  H(u)  G{uj)  =  0.  (4.1) 

If  the  family  of  bases  {ipn(x)}  is  complete  and  linearly  independent,  but  not 
orthonormal,  the  wavelet  transform  is  called  a  biorthogonal  wavelet  transform. 
Biorthogonal  wavelets  have  dual  basis  functions.  More  generally,  if  the  family  of 
wavelets  {ipn{x)}  is  not  linearly  independent  (redundant)  and  overcomplete,  they 
may  form  a  wavelet  frame  representation  [64]. 

For  a  dyadic  wavelet  transform,  the  orthonormal  constraint  is  relaxed,  so  we  can 
have  distinct  decomposition  and  reconstruction  wavelets  as  long  as  the  corresponding 
low-pass  filter  H(lj)  and  high-pass  filters  G(uj),K(uj)  satisfy  [51] 

\H{u)\2  +  G(w)K{u)  =  1.  (4.2) 

The  above  discussion  is  for  one  dimension  functions.  It  can  easily  extend  to  two 
dimensions  through  a  few  well  known  methods  [13,  51].  In  this  section,  we  focus  on 
continuous  sub-octave  wavelet  transforms.  We  first  discuss  one-dimensional  multi- 
scale  sub-octave  wavelet  transforms  and  corresponding  sub-octave  wavelet  represen- 
tations. We  then  introduce  two-dimensional  sub-octave  wavelet  transforms  (SWT) 
and  2-D  sub-octave  wavelet  representations. 

4.2.1     One-Dimensional  Sub-Octave  Wavelet  Transform 

If  we  further  divide  an  octave  frequency  band  into  M  equally-spaced  sub-octave 
bands  (here,  M  is  limited  to  a  power  of  two),  then  M  wavelet  bases  can  be  used  to 
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capture  the  detail  information  of  a  signal  in  each  sub-octave  frequency  band.  The 
M  wavelet  functions  are  represented  as  ipm(x)  G  L2(R),  where  m  =  {1,2,  ..,M}, 
L2(R)  denotes  the  space  of  measurable,  square-integrable  1-D  real  functions.  An  M 
sub-octave  wavelet  transform  of  a  1-D  function  f(x)  G  L2(R)  at  scale  2J  (level  j) 
and  location  x,  and  for  an  ra-th  sub-octave  frequency  band  is  defined  as 


/+oo 
f(t)f/%(x-t)dt,  (4.3) 

-oo 


J 


where  ip%j(x)  =  ^j4'm{^j)  is  the  dilation  of  the  m-th.  wavelet  basis  ipm(x),  at  scale  2\ 
m  =  {1, 2, ...,  A/},  and  j  G  Z.  In  the  frequency  domain,  we  can  write 

Wpf(co)  =  />)^(2M,  (4-4) 

by  taking  the  Fourier  transform  of  Equation  (4.3).  A  scaling  approximation  of  a  1-D 
function  f(x)  is  defined  as 

Sjf(x)  =  /*Mz)-  (4-5) 

To  provide  a  more  flexible  choice  for  the  M  basis  wavelets,  we  impose  that  the  wavelet 
functions  satisfy  a  wavelet  frame  condition  (similar  to  Mallat  and  Zhong  [51]) 

+oo        A/ 

^  <  E  E  l^ro(2V)l2  <  5, 

j  —  —  OO    771=1 

where  A  and  B  are  positive  constants  and  uj  G  R.  In  the  spatial  domain,  we  have 


+oo        M 

A\\m\\2  <  e  e  iiwr/wir  <  ^ii/wii2- 

j  =— oo  m=l 
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The  function  f(x)  can  be  recovered  from  its  sub-octave  wavelet  transform  by  the 
formula 

+00        M  +00        A/ 

/(*)=  E  £wt/*t5(*)  =  E  £/*^*t£(*),        (4-6) 

j=— 00  m=l  j  —  —  oo  m=l 

where  7m(a;)  is  the  m-th  synthesis  wavelet  for  the  inverse  sub-octave  wavelet  trans- 
form. The  set  of  frequency  response  of  {tp!£(x)  \m  =  1,2,...,  M }  together  at  any  scale 
2J  are  required  to  capture  the  details  within  an  octave  frequency  band.  Finally,  for 
perfect  reconstruction  of  f(x),  analysis  wavelets  ip^(x)  and  synthesis  wavelets  7^(2:) 
should  satisfy 

+00        M 

£     £^"(2'u,)7m(2M  =  1.  (4.7) 

j——oo  m=l 

Equation  (4.7)  can  be  obtained  by  taking  the  Fourier  transform  on  both  sides  of 
Equation  (4.6).  To  ensure  exact  reconstruction,  the  frequency  axis  is  covered  by 
both  analysis  and  synthesis  wavelets.  Thus,  the  wavelets  7m(x)  can  be  any  functions 
whose  Fourier  transforms  satisfy  Equation  (4.7).  There  are  certainly  many  choices 
for  analysis  and  synthesis  wavelets  that  satisfy  Equation  (4.7).  For  de-noising  and 
feature  enhancement  purposes,  we  are  interested  in  the  class  of  wavelets  which  are 
an  approximation  to  the  first  or  second  order  derivatives  of  a  smoothing  function, 
such  as  spline  functions  of  any  order.  A  1-D  sub-octave  wavelet  transform  can  be 
easily  extended  to  2-D  by  computing  sub-octave  wavelet  coefficients  along  horizontal 
and  vertical  directions  [51],  as  explained  next  in  the  following  section.  Extensions  to 
higher  dimensions  are  straight  forward  and  analogous. 

4.2.2     Two-Dimensional  Sub-Octave  Wavelet  Transform 

Daubechies  described  two  ways  to  extend  1-D  orthonormal  wavelet  transforms 
to  two  dimensions  [13].  Here  we  adopt  the  way  of  dyadic  wavelet  extension  to  two 
dimensions  introduced  by  Mallat  and  Zhong  [51]  by  computing  sub-octave  wavelet 
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coefficients  along  horizontal  and  vertical  directions.  An  M  sub-octave  wavelet  trans- 
form of  a  2-D  function  f(x,y)  G  L2(R2)  at  scale  2J  (level  j)  and  location  (x,  y),  for 
the  m-th  sub-octave  frequency  band  is  defined  as 


Wfmf(x,y)  =f*ifemfay),  (4.8) 

where  ^/"(a;,  y)  —  ^jipd'm{^j,  %■),  d  =  {1,2}  (for  horizontal  and  vertical  directions), 


m  —  {1,  2, ...,  M },  and  j  G  Z.    L2{R  )  denotes  the  space  of  measurable,  square- 
integrable  2-D  functions.  In  the  Fourier  domain,  Equation  (4.8)  simply  becomes 

Wf>mf(ux,ujy)  =  f(u;x:ujy)^m(2^x,2^y).  (4.9) 

The  function  f(x,y)  can  be  recovered  from  its  2-D  sub-octave  wavelet  transform  by 
the  formula 

+  oo         2        A/ 

/(*,»)    =      E    E  E  <""/*7jm(x,!/) 

j  =  — oo  d=l  m=l 

+oo         2        A/ 

E  E  E  /*!$m  *-#"(*,*)■  (4-io) 

J  —  —  00    d=l     771=1 

For  perfect  reconstruction,  ijj7^(x,y)  and  ^(x,y)  must  satisfy 

+oo         2        M 

E     E  E  ^(2^,2^)7^(2^,2^)  =  1.  (4.11) 

j=  —  oo    d=\   m=l 

This  exact  reconstruction  condition  is  obtained  by  taking  the  Fourier  transform  of 
Equation  (4.10). 
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4.3     Discrete  Sub-Octave  Wavelet  Transform 

Continuous  wavelet  transforms  are  useful  to  demonstrate  the  properties  of  wavelet 
decomposition  and  reconstruction  and  are  helpful  for  theoretical  approval  of  the  per- 
fect reconstruction  while  discrete  wavelet  transforms  are  practical  important  for  dis- 
crete signal  and  digital  image  processing.  The  transform  parameters  in  a  sub-octave 
wavelet  transform  are  continuous  variables.  This  results  in  a  highly  redundant  rep- 
resentation. It  is  possible  to  discretize  these  parameters  and  still  achieve  perfect 
reconstruction  [64].  For  digital  image  processing,  the  sub-octave  wavelet  transform 
of  a  discrete  function  can  be  carried  out  through  uniform  sampling  of  a  continu- 
ous sub-octave  wavelet  transform.  Below,  we  describe  the  discrete  formulation  of  a 
sub-octave  wavelet  transform. 

4.3.1     One-Dimensional  Discrete  Sub-Octave  Wavelet  Transform 

In  the  discrete  domain,  scales  are  also  discrete  and  limited  by  the  finest  scale  of 
one  unit.  A  sub-octave  wavelet  transform  can  similarly  be  decomposed  into  dyadic 
scales  and  we  can  get  a  perfect  reconstruction  through  its  corresponding  inverse 
sub-octave  wavelet  transform.  In  general,  a  function  can  be  decomposed  into  fine- 
to-coarse  (dyadic)  scales  by  its  convolution  with  dilated  wavelets  {ip2j{x)}jGz-  This 
can  be  done  through  repeated  smoothings  (low  pass  filtering)  and  detail  finding  (high 
pass  filtering).  In  the  discrete  domain,  because  of  the  limitation  of  the  finest  scale, 
scales  have  to  be  greater  than  or  equal  to  1,  so  we  let 


+oo     M 

l^)|2  =  E  E^(2M7m(2M-  (4-12) 

j=l  m=l 
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Thus,  for  a  J-level  discrete  sub-octave  wavelet  transform,  we  can  write 

J        M 

I^MI2  =  E  E  ^(2M  7m(2M  +  mJuj)\2.  (4.13) 

j=\  rn—\ 

The  notation  {Wpf(n),  Sjf(n)  \j  =  1,2, ...,  J,  and  r?i  =  1,2,  ...,M}  is  defined  as 
the  wavelet  representation  of  a  discrete  function  /(n)  under  a  J-level  discrete  M 
sub-octave  wavelet  transform.  W™  f(n)  and  Sjf(n)  are  uniform  samplings  of  their 
continuous  counterparts  respectively. 

If  we  let  J  =  1,  then  Equation  (4.13)  becomes 


M 

mco)\2  =  £  ^m(2u)  7m(2u;)  +  |^(2u;)|2.  (4.14) 

m=l 


Let  us  assume  that  for  each  basis  wavelet  pair  ipm(ij)  and  7m(u;),  there  exists  a  pair  of 
corresponding  functions  Gm(uj)  and  Km{uj)  which  need  to  be  determined  and  (with 
certain  temporal  shift  t)  satisfy 

^pm{2uj)  =  £{u)Gm(uj)e-ltuJ1 
7m(2u;)  =  fi{uj)Km{uj)eltu). 

And,  for  scaling  function  (p{u)),  there  exists  a  function  H(u>)  which  satisfies 

ip(2u)  =  <p{u)  H  (u)  e~itw. 

Replacing  xf;m(2u),  7m(2u;),  and  <p(2uj)  in  Equation  (4.14),  we  obtain  a  fundamental 
relation  for  sub-octave  wavelet  transforms  (SWT);  that  is, 


M 

|//(u;)|2  +  E  Gm{u)  Km{v)  =  1-  (4-15) 

m=l 
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Figure  4.1.     A  3-level  SWT  decomposition  and  reconstruction  diagram  for  a  1-D 
function. 

The  discrete  sub-octave  wavelet  transform  of  a  function  f(n)  G  l2(Z)  can  be 
implemented  by  using  the  following  recursive  relations  between  two  consecutive  levels 
j  and  j '  +  1 

W?+J(u)  =  Sjf(u)Gm(Vu),  (4.16) 


Sj+1f(u;)  =  Sjf(uj)H(yu), 


(4.17) 


where  j  >  0,  1  <  m  <  M,  and  S0f(uj)  —  f{oj).  And,  the  reconstruction  Sof(cu) 
from  a  sub-octave  wavelet  decomposition  can  be  implemented  through  the  recursive 
relation 

M 

(4.18) 


S3f(co)  =  Sj+1f(u)  H(Vu)  +  £  W?+1f(u)  Km(2^). 
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Figure  4.2.  Divisions  of  the  frequency  bands  under  the  SWT  shown  in  Figure  4.1. 

where  H  is  the  complex  conjugate  of  H.  A  three-level  M  sub-octave  wavelet  decom- 
position and  reconstruction  process  based  on  the  above  recursive  relations  is  shown  in 
Figure  4.1.  The  corresponding  divisions  of  frequency  bands  are  depicted  graphically 
in  Figure  4.2.  In  general,  for  an  M  sub-octave  analysis  and  synthesis,  we  require 
M  pairs  of  corresponding  basis  wavelets.  A  SWT  is  a  multiwavelet  transform  with 
a  single  scaling  function  [60,  61].  When  M  is  a  power  of  2,  we  can  carry  out  de- 
composition and  reconstruction  using  a  set  of  FIR  filters  corresponding  to  a  single 
pair  of  basis  wavelets.  Figure  4.3  presents  a  filter  bank  for  carrying  out  a  2-level 
4  sub-octave  decomposition  and  reconstruction  using  FIR  filters  corresponding  to  a 
single  pair  or  two  pairs  of  basis  wavelets.  This  describes  a  more  general  way  to  do  the 
sub-octave  decomposition  where  fine-to-coarse  octave  decomposition  and  sub-octave 
band  splitting  can  be  carried  out  through  two  sets  of  different  FIR  filters.  It  reduces 
to  the  case  [42]  when  H  =  Hs  and  G  =  Gs  where  Hs  and  Gs  are  used  for  sub-octave 
band  splitting. 

4.3.2     Two-Dimcnsional  Discrete  Sub-Octave  Wavelet  Transform 

For  the  decomposition  of  a  2-D  discrete  function,  we  let  the  frequency  response 
of  a  scaling  function  be  defined  in  the  formula 


+oo     2        M 

\<p(ujx,uy)\2  =  E  E  E  V/i/m(2Vc,  Vuy)  7*m(2^,  2Ju;y).  (4.19) 

j=l  d=\  m=\ 
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Figure  4.3.  A  2-level  4-sub-octave  decomposition  and  reconstruction  of  a  SWT. 
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Figure  4.4.  Frequency  plane  tessellation  and  filter  bank,  (a)  Division  of  the  frequency 
plane  for  a  2-level,  2-sub-octave  analysis,  (b)  Its  filter  bank  along  the  horizontal 
direction. 
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For  a  J-level  2-D  discrete  sub-octave  wavelet  transform,  we  can  formulate 

J       2       M 

|£K,u;y)|2  =  E  E  E  ^d'm(2j^,  2Ju,y)  7d'm(2-^,  Vuy)  +  l^w,,  2J^)|2. 

j|=l  d=l  m=l 

(4.20) 
If  the  scaling  approximation  of  a  function  f(x,y)  at  scale  2J'  is  represented  by 

Sjf(x,y)  =  f*<p3f(xty),  (4.21) 

then  {  Wf'm  f{nx,ny),Sjf{nx,ny)  \d  =  1,2,  j  =  1, ..,  J,  and  ra  =  1,..,A/}  is  called 
the  wavelet  representation  of  a  discrete  function  f{nx,ny)  for  a  2-D  J-level  discrete 
M  sub-octave  analysis.  In  general,  <p(x,y)  is  a  2-D  scaling  function  and  ipd,m(x,y) 
and  jd,m(x,y)  are  the  ra-th  directional  analysis  and  synthesis  wavelets.  There  are 
many  choices  of  these  functions  that  satisfy  Equation  (4.20).  Similar  to  the  way  2-D 
wavelets  are  constructed  using  1-D  wavelets  [51],  we  use  tensor  products  to  construct 
2-D  sub-octave  wavelets  using  1-D  sub-octave  wavelets.  Thus  we  can  write 

ft>m(Vujx,  Vuiy)  =  $m{yu)x)w(y-lu)y),  (4.22) 

j>2>m(2jux,  2juy)  =  <p{2j-1ux)i)m(2'ujy),  (4.23) 

ip(2jux,  2juy)  =  0{2jujx)<p{2juy).  (4.24) 

Through  this  construction,  we  implemented  a  2-D  sub-octave  wavelet  transform  using 
1-D  convolution  with  FIR  filters  of  the  1-D  sub-octave  wavelet  transform  previously 
described.  Figure  4.4(a)  shows  the  division  of  the  frequency  plane  under  a  2-level 
SWT  where  M  =  2.  Figure  4.4(b)  shows  the  corresponding  filter  bank  along  the 
horizontal  direction  where  the  curve  shown  in  red  corresponds  to  the  analytic  fil- 
ter ip1,l(2ux,2ujy)  (Wi'  )  projected  along  the  lox  axis,  the  black  curve  for  W\" .  the 
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Figure  4.5.  Smoothing,  scaling,  and  wavelet  functions  for  a  SWT.  These  functions 
are  used  for  a  2-sub-octave  analysis. 

magenta  curve  for  W2'  ,  the  blue  curve  for  W2  ,  and  the  green  one  for  52-  A  2-D 
sub-octave  wavelet  transform  can  be  implemented  by  1-D  convolution  with  FIR  filters 
along  horizontal  and  vertical  directions.  The  details  along  the  diagonal  directions  are 
embedded  in  the  details  along  horizontal  and  vertical  directions.  In  Figure  4.5,  a 
fourth-order  spline  smoothing  function,  its  first  and  second  derivative  approximation 
as  sub-octave  wavelets,  and  the  corresponding  scaling  function  are  shown. 

A  dyadic  wavelet  transform  can  be  a  special  case  of  a  sub-octave  wavelet  trans- 
form with  M  =  1.  As  an  example,  a  discrete  2-sub-octave  wavelet  transform  is  shown 
to  divide  the  details  of  an  octave  band  into  details  of  2  sub-octave  bands.  As  shown 
in  Figure  4.6,  one  sub-octave  can  be  used  for  detecting  local  maxima  while  the  other 
sub-octave  band  identifies  zero-crossings.  The  dashed  curve  corresponds  to  the  fre- 
quency response  of  a  first  order  derivative  approximation  of  a  smoothing  function 
and  the  dash-dot  curve  shows  the  frequency  response  of  a  second  order  derivative 
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Figure  4.6.  An  example  of  level  one  analytic  filters  for  2  sub-octave  bands  and  a 
low-pass  band.  The  dashed  curve  is  the  frequency  response  of  a  first  order  deriva- 
tive approximation  of  a  smoothing  function  and  the  dash-dot  curve  is  the  frequency 
response  of  a  second  order  derivative  approximation.  The  solid  curve  is  a  scaling 
approximation  at  level  one. 

approximation.  The  solid  curve  is  a  scaling  approximation  for  analysis  at  level  one. 
Thus,  these  analysis  filters  take  advantage  of  both  local  maxima  and  zero-crossing 
representations  to  characterize  local  features  emergent  within  each  scale. 

4.4     SWT-Based  De-Noising 

Our  experiments  showed  that  a  S WT  with  first  and  second  order  derivative  ap- 
proximation of  a  smoothing  function  as  its  basis  wavelets,  separated  coefficients  best 
characterized  by  feature  energy  from  coefficients  characterized  by  noise  energy. 

Sub-octave  wavelet  coefficients  can  be  modified  by  hard  thresholding  for  noise 
reduction  by 

W/'m'7(x)  =  Mdfl  TH{Wf>mf(x)/M$m,  t*m),  (4.25) 


fd,m 


d,rn 


M$m  =  max(\W°>mf(x)\), 


(4.26) 
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where  d  =  {1,2}  (omitted  for  1-D  signals),  j  =  {1, ...,  k},  and  k  <  J,  m  =  {1, ...,  Af }, 
and  £2]m  is  (in  general)  a  threshold  related  to  noise  and  scale.  The  processed  result 
Wj'Ti'*f(x)  is  a  modified  coefficient.  Position  x  denotes  n  for  1-D  signal  processing 
and  (nx,ny)  for  2-D  image  noise  reduction. 

SWT  coefficients  can  be  processed  through  soft  thresholding  wavelet  shrinkage 
for  noise  reduction  by 

W/'m'7(x)  =  M%m  Ts{W*'mf{x)/M%m,  4"),  (4-27) 

M^  =  max(|W/'m/(x)|),  (4.28) 

where  d  =  {1,  2}  (omitted  for  1-D  signals),  j  =  {1, ...,  h},  k  <  J,  m  =  {1, ...,  M},  and 
t2jm  is  a  noise  and  scale-related  threshold.  Again,  the  result  W,  'm'*/(x)  is  a  processed 
coefficient.  Similarly,  position  x  denotes  n  for  1-D  signal  processing  and  (nx,ny)  for 
2-D  image  noise  reduction.  Recall  that  Donoho  and  Johnstone's  soft  thresholding 
method  used  a  single  global  threshold  [19].  However,  since  noise  coefficients  under 
a  SWT  have  average  decay  through  fine-to-coarse  scales,  we  regulate  soft  threshold- 
ing wavelet  shrinkage  by  applying  coefficient-dependent  thresholds  decreasing  across 
scales  [72].  When  features  and  noise  can  be  clearly  separated  at  the  finer  levels  of 
scale,  applying  hard  thresholding  instead  of  soft  thresholding  may  further  improve 
performance.  However,  hard  thresholding  may  fail  to  smooth  noise  if  the  noise  is 
strong  locally. 
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Table  4.1.     Quantitative  measurements  of  performance  for  de-noising  and  feature 
restoration. 


Method  or  Measurement 

Blocks 

Bumps 

HeaviSine 

Doppler 

Noisy  Signal:  ojo^ 

6.856 

6.735 

6.895 

7.017 

Noisy  Signal:  10 logl0{a2s/tf)     (dB) 

16.726 

16.567 

16.771 

16.923 

Restored  Result:  osjan 

27.258 

20.988 

35.453 

20.129 

Restored  Result:  10logl0(a2Ja2)  (dB) 

28.779 

26.439 

30.993 

26.076 

4.5     SWT-Based  Enhancement  with  Noise  Suppression 

Through  a  nonlinear  enhancement  function  or  a  generalized  adaptive  gain  opera- 
tor, SWT  coefficients  were  modified  for  contrast  enhancement  by 


rd,m,*  fl  \   __    -\/rd,m   j-.  f    ™  j       J\nxiny) 


W;^f{nxlny)=Ma2rEOP 


Md.m 
2i 


rd,m 


d.m 


M%m  =  mzK(\W°>mf(nx,ny)\) 


(4.29) 


(4.30) 


where  d  =  {1,  2},  1  <  m  <  M,  and  1  <  j  <  J.  The  pointwise  operator's  output, 
Wj'm'*f(nx,ny)  is  simply  a  processed  coefficient.  The  enhancement  operator  Eop 
can  be  ENL  or  Eqag  presented  in  Section  2.4.  Since  these  operators  are  defined 
on  input  values  between  -1  to  1,  we  normalize  sub-octave  wavelet  coefficients  before 
applying  the  operators. 

4.6     Application  Samples  and  Analysis 

In  this  section,  we  present  several  experimental  results  of  SWT  based  de-noising 
and  enhancement  for  both  1-D  signal  and  2-D  images.  First  we  present  experimen- 
tal results  based  on  a  1-D  sub-octave  wavelet  transform.  We  show  the  de-noising 
capability  of  our  method  with  feature  restoration  compared  to  existing  de-noising 
methods. 
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Figure  4.7.  De-noised  and  restored  features  from  the  SWT-based  algorithm.  From 
top  to  bottom:  original  signal;  noisy  signal;  de-noised  signal;  overlay  of  original  and 
de-noised  signal. 
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Experimental  results  of  SWT-based  de-noising  are  shown  in  Figure  4.7.  In  the 
fourth  plot  shown  in  each  Figure  4. 7(a)- (d),  the  de-noised  results  are  overlaid  with 
their  corresponding  original  signals.  Table  4.1  shows  quantitative  measurements  of 
each  result  shown  in  Figures  4.7.  In  comparison  to  previously  published  methods 
processing  the  same  signals  [8],  results  of  the  SWT-based  method  are  noticeably  im- 
proved and  basically  free  from  artifacts,  pseudo-Gibbs  phenomena.  In  the  next  exper- 
iment, we  show  the  limitations  of  a  traditional  D WT  for  characterizing  band-limited 
high  frequency  features.  Figures  4.8(a)-(b)  show  the  original  and  noisy  "Doppler" 
signals  with  their  corresponding  5-level  D WT  and  a  coarse  scale  approximation.  The 
finest  scale  band-limited  features  (see  the  second  plot  in  Figure  4.8(a))  are  weak  and 
obscured  when  noise  is  present  (see  the  second  plot  in  Figure  4.8(b)).  These  high  fre- 
quency band-limited  features  are  lost  in  a  soft  thresholding  de-noising  method,  shown 
in  Figure  4.9(b).  Figures  4.8(c)-(f)  show  2-sub-octave  wavelet  transforms  of  the  orig- 
inal and  noisy  "Doppler"  signals.  Figures  4.8(c)  and  4.8(e)  show  first  sub-octave 
coefficients  while  Figures  4.8(d)  and  4.8(f)  display  second  sub-octave  coefficients. 
Figure  4.9  shows  de-noised  and  enhanced  results  of  noisy  "Doppler"  under  a  DWT 
and  a  SWT.  The  loss  of  high  frequency  band-limited  features,  made  the  result  from 
the  DWT-based  method  less  attractive  than  the  SWT-based  method,  for  processing 
medical  images,  such  as  microcalcification  in  mammograms. 

Next,  experimental  results  of  enhancement  with  noise  suppression  of  mammo- 
graphic  images  via  a  2-D  sub-octave  wavelet  transform  are  presented.  We  demon- 
strate the  advantage  of  enhancement  without  amplifying  noise,  including  background 
noise  [37].  In  the  first  experiment  for  image  enhancement,  we  attempt  to  enhance 
the  visibility  of  a  radiographic  image  of  a  RMI  (Radiation  Measurements  Inc.,  Mid- 
dleton,  WI)  model  156  phantom  with  simulated  masses  embedded.  Figure  4.10(a) 
presents  a  low  contrast  image  of  the  RMI  model  156  phantom  with  simulated  masses 
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Figure  4.8.    Limitations  of  a  DWT  for  characterizing  band-limited  high  frequency 
features. 
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Figure  4.9.  De-noised  and  enhanced  results  of  a  noisy  "Doppler"  signal  under  a  DWT 
(25.529dB)  and  a  SWT  (26.076dB). 
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(b) 

Figure  4.10.  Enhancement  with  noise  suppression,  (a)  A  low  contrast  image  of  RMI 
model  156  phantom  with  simulated  masses  embedded,  (b)  Enhancement  by  tradi- 
tional histogram  equalization,  (c)  SWT-based  enhancement  with  noise  suppression. 

of  distinct  sizes.  Figure  4.10(b)  shows  enhancement  by  traditional  histogram  equal- 
ization. SWT-based  enhancement  with  noise  suppression  is  shown  in  Figure  4.10(c). 
Unsharp  masking,  a  popular  technique  for  enhancing  radiographs,  failed  to  enhance 
barely  seen  masses  (low  frequency  features)  in  Figure  4.10(a).  In  comparison  to  tra- 
ditional histogram  equalization,  the  SWT-based  method  enhanced  the  visibility  of 
masses  without  amplifying  noise. 

In  the  next  experiment,  we  enhance  low  contrast  mammographic  images  contain- 
ing a  microcalcification  cluster.  Figure  4.11(a)  shows  a  region  from  a  low  contrast 
mammogram  containing  a  distinct  microcalcification  cluster.  Next,  enhancement  by 
traditional  unsharp  masking  is  presented  in  Figure  4.11(b).  Figure  4.11(c)  shows  the 
result  of  SWT-based  enhancement  with  noise  suppression.  In  practice,  a  radiologist 
may  want  to  view  certain  suspicious  areas  of  low  contrast  within  a  large  digital  mam- 
mogram for  potential  breast  lesions  with  a  magnifier  to  improve  visibility  of  an  area. 
In  Figure  4.11(d),  we  try  to  provide  a  similar  function  by  improving  the  visibility  of  a 
local  region  of  interest  (ROI)  through  SWT-based  enhancement  with  noise  suppres- 
sion. The  region  within  the  black  square  (120x120)  is  enhanced.  Note  that  the  area 
does  not  have  to  be  a  square,  but  a  rectangle.   Traditional  unsharp  masking  shown 
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in  Figure  4.11(b)  enhances  the  area  of  the  microcalcification  cluster  slightly  but  also 
amplifies  noise.  As  shown  in  Figures  4.11(c)  and  4.11(d),  enhancement  under  a  SWT 
makes  barely  seen  microcalcification  clusters  more  visible  without  amplifying  noise. 
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(a) 


(b) 


(c) 


(d) 


Figure  4.11.  Enhancement  with  noise  suppression,  (a)  Area  of  a  low  contrast  mam- 
mogram with  a  microcalcification  cluster,  (b)  Best  enhancement  by  traditional  un- 
sharp  masking,  (c)  SWT-based  enhancement  with  noise  suppression,  (d)  SWT-based 
enhancement  of  a  local  region  of  interest  (ROI)  with  noise  suppression. 


CHAPTER  5 
PERFORMANCE  MEASUREMENT  AND  COMPARATIVE  STUDY 

In  this  chapter,  we  present  several  measures  of  relevant  quantitative  metric  for 
evaluating  an  algorithm's  performance  and  show  a  few  comparative  studies  of  quan- 
titative measurements  between  our  algorithms  and  other  published  methods.  In 
Section  5.1,  a  few  quantitative  measurements  used  in  this  dissertation  are  described 
and  formulated  mathematically.  Section  5.2  shows  the  quantitative  results  of  our 
algorithms  and  other  relevant  methods  for  signal/image  restoration.  In  Section  5.3, 
we  present  quantitative  measurements  of  image  enhancement  among  our  developed 
algorithms  and  a  few  other  related  methods.  Earlier  chapters  demonstrated  visual 
quality  of  de-noised  as  well  as  enhanced  results  (visual  performance)  while  this  chap- 
ter focuses  on  quantitative  performance. 

5.1     Performance  Metric  for  Quantitative  Measurements 

The  quality  of  a  noisy  signal/image  is  often  measured  by  the  ratio  of  signal  vari- 
ance to  noise  variance  (or  signal  energy  level  of  variation  to  noise  energy  level)  using 
a  log  scale.  The  quantity  is  called  signal-to-noise  ratio.  A  signal/image  is  more  likely 
severely  degraded  when  noise  level  is  high  (low  signal-to-noise  ratio).  We  have  used 
the  quantitative  term  signal-to-noise  ratio  when  displaying  our  earlier  experimental 
results.  A  formal  definition  of  signal-to-noise  ratio  is  given  below. 

Signal-to-noise  ratio  (SNR)  is  defined  as 


2 

SNR  =  10  lognfy, 
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where  a2  and  a2  are  signal  variance  and  the  average  noise  energy  (average  least 
squares  error  between  a  signal  and  its  noise-free  original  version),  respectively  [30]. 
Here,  we  denote  an  ideal  (original)  signal  as  g  and  the  restored  version  from  its  noise 
signal  as  g*.  For  1-D  signals,  a2  is  defined  as 


iy    n=l 


and  a2  is  defined  by 


1     N 
9=M  S#(n), 

iV      71=1 


tf  =  ^lL(9*(n)-g(n)y 

7V      71=1 


For  2-D  images,  a2  is  defined  as 


M       N 

IV1    iV      771=1    71=1 


M       N 

2  =  ^77  £  E9(m,n) 

IVI    JV       771=1    71=1 


while  a2  is  defined  by 


i  M      N 

a^  =  WN   E   £(#*(m>n)  -g{m,n))' 

1V1    iV      771=1    71=1 


The  term  SNR  was  frequently  used  in  many  noise  reduction  algorithms  for  quan- 
titative measurements  of  performance  and  has  also  been  applied  to  measure  our 
de-noising  methods.  The  quality  improvement  of  a  signal/image  can  be  measured  by 
an  improved  (higher)  SNR  compared  to  the  SNR  of  its  noisy  version. 

The  quantitative  measurements  of  average  square  errors  were  used  to  measure  the 
performance  of  noise  reduction  algorithms  [19,  20,  8].  For  our  comparative  studies,  we 
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use  the  same  measurements  to  compare  the  de-noising  performance  between  our  algo- 
rithms and  other  published  methods,  such  as  Donoho  and  Johnstone's  thresholding- 
based  wavelet  shrinkage,  Coifman  and  Majid's  wavelet  packet-based  de-noising  [9], 
and  Chen  and  Donoho's  basis  pursuit  de-noising  [5].  Average  square  error  (ASE)  for 
1-D  signal  recovery  is  defined  as 


1 1   *         ii- 

ACl-i         \\9      ~  9\\2,N 

S         N ' 


where  H/H^w  =  X^=i(/(n))2-  Average  square  error  (ASE)  for  2-D  image  restoration 

is  defined  by 

1 1 /i*  -  nil2 
AQj?  _    I  iff  9\\2,{M,N} 

ASE  -  MN ' 

where  \\f\\l{MjN)  =  E^E^/Kn))2. 

Another  often  used  metric  for  measuring  an  algorithm's  performance  is  the  root 
mean  squared  error  [8].  We  also  use  the  quantity  for  measuring  our  algorithm's 
performance  compared  to  Coifman  and  Donoho's  methods.  The  root  mean  squared 
error  (RAISE)  for  1-D  signal  restoration  can  be  defined  as 


RMSE  =  ^,(r(n)-g(n)P 


For  2-D  image  processing,  the  root  mean  squared  error  (RMSE)  for  image  restoration 
is  defined  by 

RMSE  =  i/^™=1  ^=1  fo*(m'  n)  ~  9^m'  rc))2 


/  MN 

For  image  enhancement,  contrast  has  been  used  as  a  major  measurement  for 
an  algorithm's  performance.  In  order  to  perform  quantitative  comparison,  we  have 
applied  the  same  quantitative  measures,  including  a  version  of  the  optical  definition 
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of  contrast  introduced  by  Morrow  et  al.  [53]  and  a  contrast  improvement  index  [39]. 
The  contrast  C  of  an  object  of  interest  is  defined  by 


C  = 


if  -  ih 


if  +  h 


where  if  is  the  mean  gray-level  (intensity)  value  of  a  particular  object  in  a  region  of 
interest,  called  the  foreground,  and  %  is  the  mean  intensity  value  of  a  surrounding 
region  called  the  background.  The  definition  of  C  above  takes  the  absolute  value 
of  the  definition  by  Laine  et  al.  [39].  The  reason  is  that  sometimes  an  object  of 
interest  may  have  lower  average  intensity  than  its  surrounding  area  and  the  techniques 
introduced  here  can  be  applied  to  other  image  processing  applications.  We  avoid 
having  a  negative  value  for  the  contrast  of  an  object  by  using  this  expression.  A 
contrast  improvement  index  (CII)  can  be  used  as  a  quantitative  measure  of  contrast 
improvement  and  is  defined  by 

pijj ^Processed 

^n     —     7^  > 

^Original 

where  CproCessed  and  Coriginai  are  the  contrasts  for  a  region  of  interest  in  the  enhanced 
and  original  images,  respectively.  We  use  contrast  and  CII  for  comparing  the  perfor- 
mance of  a  few  image  enhancement  methods  in  Section  5.3. 

The  performance  of  related  algorithms,  such  as  signal/image  restoration  and  im- 
age enhancement,  can  be  relatively  compared  using  similar  quantitative  measure- 
ments. The  quantities  allow  comparison  of  the  relative  performance  of  the  algorithms. 
In  the  next  two  sections,  we  will  evaluate  the  performance  of  several  related  methods 
for  noise  reduction  and  image  enhancement. 
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Table  5.1.  Quantitative  measurements:  Average  Square  Errors  -  for  various 

signal  restoration  methods. 


TV 


Method  or  Measurement 

Blocks 

Bumps 

HeaviSine 

Doppler 

WIWIn/n 

81.211 

57.665 

52.893 

50.348 

with  noise 

1.047 

0.937 

1.008 

0.999 

ideal  wavelets 

0.097 

0.111 

0.028 

0.042 

threshold  X*N 

0.395 

0.496 

0.059 

0.152 

threshold  (2  log  TV)1/2 

0.874 

1.058 

0.076 

0.324 

WP  coherent  de-noising 

0.433 

0.532 

0.133 

0.257 

xwpl  and  WB 

0.081 

0.182 

0.056 

0.126 

xwpl  and  BOB 

0.081 

0.239 

0.084 

0.104 

BPDN-JumpWave 

0.394 

0.761 

0.072 

0.282 

BPDN-HS 

0.060 

failed 

0.100 

0.319 

DWT  soft  thresholding 

0.149 

0.274 

0.049 

0.209 

D WT  de-noising  and  enhancement 

0.039 

0.138 

0.042 

0.137 

S WT  de-noising  and  enhancement 

0.065 

0.111 

0.038 

0.120 

5.2     Quantitative  Comparison  of  Signal/Image  Restoration 

This  section  presents  quantitative  measurements  for  comparing  the  performance 
of  several  algorithms  aimed  at  data  restoration.  In  Chapter  3,  we  have  used  SNR 
to  show  an  algorithm's  quantitative  performance  in  a  few  experiments  while  showing 
its  visual  performance,  such  as  in  Figures  3.11  and  3.12.  Table  4.1  in  Chapter  4 
has  also  shown  SNR  for  noisy  signals  and  restored  results  after  noise  reduction 
and  feature  enhancement.  Table  5.1  shows  the  quantitative  measurements  of  our 
algorithm  performance  compared  to  Donoho  and  Johnstone's  threshold-based  de- 
noising  methods,  Coifman  and  Majid's  wavelet  packet  based  methods,  as  well  as 
Chen  and  Donoho's  basis  pursuit  based  methods. 

The  quantity  H/H^/v/^  in  Table  5.1  is  the  average  signal  energy  per  pixel,  where 
1 1 /1 1 2,7V  =  Sn=i  f(n)2-  The  second  row  of  the  table  shows  average  square  error  of  the 
noisy  signal  to  its  noise-free  original.  The  results  in  the  third  row  were  produced  by 
ideal  wavelet  coefficient  selection  using  the  position  information  of  noise-free  wavelet 
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coefficients  (in  Donoho  and  Johnstone's  term,  with  the  help  of  an  oracle,  which  was 
called  "ideal  wavelets").  The  method  is  basically  using  wavelet  coefficients  of  the 
noise-free  signal  to  select  and  scale  the  corresponding  noisy  wavelet  coefficients  based 
on  the  noise  level  in  order  to  remove  noise  components.  The  next  two  rows  are 
the  results  of  Donoho's  wavelet  coefficient  thresholding  methods  under  two  different 
thresholds  [19].  The  row  under  "WP  coherent  de-noising"  shows  the  results  of  Coif- 
man's  adapted  waveform  analysis  and  de-noising  [9].  The  results  of  "xwpl  and  WB" 
and  "xwpl  and  BOB"  were  produced  by  "xwpl"  with  "Wavelet  Basis"  and  "Best 
Of  Bases"  algorithms,  respectively.  The  results  of  "BPDN-Jump Wave"  and  "BPDN- 
HS"  are  from  Chen  and  Donoho's  paper  [5]  (Figures  5.4  and  5.5)  with  the  dictionaries 
of  "Jump+Wave"  and  "HeaviSide",  respectively.  The  results  of  "DWT  soft  thresh- 
olding" were  obtained  through  multiscale  DWT  and  soft  thresholding  methods.  The 
results  of  "DWT  de-noising  and  enhancement"  are  our  earlier  results  via  wavelet 
shrinkage  and  feature  emphasis  of  DWT  coefficients.  The  last  row  shows  the  results 
from  our  sub-octave  wavelet  based  algorithm  for  de-noising  and  enhancement.  In 
the  table,  we  include  the  results  from  other  published  methods  only  for  comparison 
purposes. 

Table  5.2  shows  the  quantitative  measurements  of  performance  among  our  devel- 
oped algorithms,  Coifman  and  Donoho's  "Translation-Invariant  De-Noising"  methods 
[8],  and  Coifman's  wavelet  packet  based  methods  [9].  The  first  seven  rows  in  Table 
5.2  are  the  results  of  Coifman  and  Donoho's  thresholding  and  different  translation- 
invariant  de-noising  methods.  The  next  three  rows  are  from  Coifman's  wavelet  packet 
based  methods  similar  to  Table  5.1,  but  in  a  different  measure.  The  next  two  rows 
are  our  earlier  DWT-based  results  which  are  the  same  as  Table  5.1  in  a  different 
quantitative  metric.  The  last  row  shows  the  results  from  our  sub-octave  wavelet 
based  algorithm  for  de-noising  and  enhancement.  These  results  are  measured  based 
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Table  5.2.  Quantitative  measurements:  RMSE  for  various  de-noising  methods. 


Method 

Blocks 

Bumps 

HeaviSine 

Doppler 

Traditional,  S8 

0.9433 

1.0040 

0.2898 

0.5352 

Cycle  Spine,  S8 

0.8901 

0.8909 

0.2830 

0.4745 

Fully  TI,  S8 

0.8458 

0.8733 

0.2855 

0.4554 

Soft,  Haar 

0.6574 

1.2291 

0.6086 

1.0881 

Soft,  Haar  TI 

0.4817 

0.8703 

0.2247 

0.7339 

Hard,  Haar 

0.2835 

0.7038 

0.4312 

0.7383 

Hard,  Haar  TI 

0.1708 

0.3966 

0.1818 

0.3893 

WP  coherent  de-noising 

0.6584 

0.7299 

0.3660 

0.5070 

xwpl  and  WB 

0.2847 

0.4273 

0.2378 

0.3563 

xwpl  and  BOB 

0.2847 

0.4895 

0.2903 

0.3237 

DWT  soft  thresholding 

0.3319 

0.5236 

0.2215 

0.4521 

DWT  de-noising  and  enhancement 

0.1819 

0.3728 

0.2050 

0.3696 

SWT  de-noising  and  enhancement 

0.2567 

0.3334 

0.1974 

0.3477 

on  "Root  Mean  Squared  Errors"  (RMSE).  Coifman  and  Donoho  [8]  claimed  they 
were  using  RMSE,  but  their  results  were  computed  using  the  root  of  overall  square 
errors  instead  (probably  caused  by  a  programming  error).  Thus,  the  difference  be- 
tween the  first  seven  rows  in  Table  5.2  and  results  in  Coifman  and  Donoho  [8]  is  that 
each  quantitative  measure  here  is  equal  to  the  corresponding  result  in  Coifman  and 


Donoho  [8]  divided  by  ^2048. 

5.3     Quantitative  Comparison  of  Image  Enhancement 

In  this  section,  we  compare  the  quantitative  performance  of  image  enhancement 
of  a  few  related  methods.  Figures  4.10  and  4.11  have  shown  the  visual  performance 
of  image  enhancement  with  noise  suppression  under  sub-octave  wavelet  transforms. 
In  the  next  few  comparative  studies,  image  enhancement  is  carried  out  on  low  con- 
trast mammographic  images  containing  a  mass,  or  simulated  mammographic  fea- 
tures. Contrast  enhancement  with  noise  suppression  for  mammographic  images  with 
or  without  additive  noise  is  shown  in  Figure  5.1.  Figure  5.1(a)  shows  an  area  of  a  low 
contrast  mammogram  with  a  lobular  mass.   Because  the  contrast  is  poor,  the  mass 
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J       L 


(a) 


(b) 


(c) 


(d) 


(e) 


(f) 


(g) 


(h) 


(0 

Figure  5.1.  Enhancement  results,  (a)  Area  of  a  low  contrast  mammogram  with  a 
mass,  (b)  Enhancement  of  (a)  by  traditional  contrast  stretching,  (c)  Enhancement  of 
(a)  by  traditional  unsharp  masking,  (d)  SWT-based  enhancement  of  (a)  with  noise 
suppression,  (e)  The  same  area  of  a  low  contrast  mammogram  contaminated  with 
additive  Gaussian  noise,  (f)  Enhancement  of  (e)  by  traditional  contrast  stretching, 
(g)  Enhancement  of  (e)  by  traditional  unsharp  masking,  (h)  SWT-based  enhance- 
ment with  noise  suppression,  (i)  Hand-segmented  mass  and  ROI  for  quantitative 
measurements  of  performance. 
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Table  5.3.  Comparison  of  contrast  values  obtained  by  traditional  contrast  stretching 
(CST),  unsharp  masking  (UNS),  and  multiscale  nonlinear  processing  of  sub-octave 
wavelet  transform  (SWT)  coefficients  of  a  mammogram  containing  a  mass  lesion. 


Contrast  values  for  CST,  UNS 
and  SWT  enhancement  techniques 

Method 

ORI 

CST 

UNS 

SWT 

First  row 

0.01920 

0.09736 

0.08390 

0.11954 

Second  row 

0.01918 

0.09295 

0.07809 

0.12802 

is  hardly  visible.  Figure  5.1(b)  shows  image  enhancement  by  traditional  contrast 
stretching  (scaling  the  low  contrast  image  linearly  to  the  full  display  range  of  grey 
levels).  Figure  5.1(c)  shows  enhancement  of  Figure  5.1(a)  through  unsharp  mask- 
ing. Figure  5.1(d)  shows  the  enhanced  result  of  Figure  5.1(a)  through  SWT-based 
enhancement  with  noise  suppression.  To  demonstrate  the  capability  of  the  method 
for  contrast  enhancement  with  noise  suppression,  we  added  Gaussian  noise  (14dB) 
to  Figure  5.1(a)  and  generated  a  noisy  image  shown  in  Figure  5.1(e).  The  additive 
noise  becomes  more  visible  in  its  linearly  enhanced  image,  through  traditional  con- 
trast stretching,  shown  in  Figure  5.1(f).  Figure  5.1(g)  shows  the  enhanced  image 
through  unsharp  masking.  Unsharp  masking  enhanced  the  mass  and  also  amplified 
the  noise.  Figure  5.1(h)  shows  the  enhanced  image  through  SWT-based  enhancement 
with  noise  suppression  and  the  result  is  pretty  close  to  the  enhanced  image  in  Figure 
5.1(d).  To  estimate  the  performance  of  various  enhancement  techniques  quantita- 
tively, we  manually  segmented  the  mass  and  ROI  (the  area  with  the  black  square) 
shown  in  Figure  5.1  (i).  The  quantitative  measurements  of  performance  for  enhance- 
ment techniques  presented  in  Figure  5.1  are  shown  in  Tables  5.3  and  5.4  based  on 
contrast  values  and  contrast  improvement  indices,  which  were  defined  and  explained 
in  Section  5.1. 
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Table  5.4.  Contrast  improvement  index  (CII)  for  enhancement  by  traditional  contrast 
stretching  (CST),  unsharp  masking  (UNS),  and  multiscale  nonlinear  processing  of 
sub-octave  wavelet  transform  (SWT)  coefficients  of  a  mammogram  with  a  mass. 


Contrast  Improvement  Index  for  CST, 
UNS,  and  SWT  enhancement  techniques 

Method 

ORI 

CST 

UNS 

SWT 

First  row 

1.00000 

5.06997 

4.36908 

6.22501 

Second  row 

0.99863 

4.84037 

4.06680 

6.66665 

Mathematical  models  of  phantoms  were  previously  constructed  in  Laine  et  al. 
[39].  Here,  we  adopt  the  same  models,  shown  in  Figure  5.2(a),  to  validate  these  en- 
hancement techniques  quantitatively  and  to  compare  the  sub-octave  based  method  to 
previous  dyadic  wavelet  processing  methods  for  improving  contrast  [39].  The  models 
included  mammographic  features,  such  as  microcalcifications,  cylindrical  and  spicu- 
lar  findings,  and  lobular  masses.  To  directly  measure  the  difference  of  performance 
between  our  enhancement  techniques  and  the  methods  reported  in  Laine  et  al.  [39], 
we  used  the  same  blended  image  shown  in  Figure  5.2(b).  For  more  details  about  how 
this  image  was  constructed,  we  refer  the  reader  to  Laine  et  al.  [39]. 

Enhancement  results  for  a  blended  image  with  simulated  mammographic  features 
are  shown  in  Figure  5.2.  Figure  5.2(a)  shows  a  phantom  image  with  the  mammo- 
graphic features.  These  features  are  blended  with  a  mammogram  M56,  as  shown  in 
Figure  5.2(b).  For  comparison  purposes,  Figure  5.2(c)  shows  the  result  after  process- 
ing the  blended  mammogram  in  Figure  5.2(b)  under  a  DWT-based  method  reported 
in  Laine  et  al.  [39].  Figure  5.2(d)  is  the  enhanced  result  produced  via  the  developed 
S WT-based  enhancement  techniques  with  noise  suppression.  Tables  5.5  and  5.6  show 
the  contrast  values  and  CII  of  the  mammographic  features  displayed  in  Figure  5.2(a). 
From  both  visual  appearance  and  quantitative  measurements,  we  observe  that  the 
result  (see  Figure  5.2(d))  obtained  via  the  SWT-based  method  is  overall  superior 
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(a) 


(b) 


(c) 


(d) 


Figure  5.2.  Phantom  enhancement  results,  (a)  Phantom  image,  (b)  Mammogram 
M56  with  blended  phantom  features,  (c)  Nonlinear  enhancement  under  a  DWT.  (d) 
SWT-based  enhancement  with  noise  suppression. 
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Table  5.5.  Comparison  of  contrast  values  obtained  by  multiscale  adaptive  gain  pro- 
cessing of  dyadic  wavelet  transform  (DWT)  and  sub-octave  wavelet  transform  (SWT) 
coefficients.  Mammographic  features:  minute  microcalcification  cluster  (MMC),  mi- 
crocalcification  cluster  (MC),  spicular  lesion  (SL),  circular  (arterial)  calcification 
(CC),  and  well-circumscribed  mass  (WCM). 


Contrast  values  for 
DWT  and  SWT  enhancement  techniques 

Feature 

MMC 

MC 

SL 

CC 

WCM 

Original 

0.050 

0.033 

0.024 

0.037 

0.012 

DWT 

0.198 

0.204 

0.165 

0.214 

0.157 

SWT 

0.342 

0.253 

0.240 

0.292 

0.112 

Table  5.6.  CII  for  enhancement  by  multiscale  adaptive  gain  processing  of  dyadic 
wavelet  transform  (DWT)  and  sub-octave  wavelet  transform  (SWT)  coefficients. 
Mammographic  features:  minute  microcalcification  cluster  (MMC),  microcalcifica- 
tion cluster  (MC),  spicular  lesion  (SL),  circular  (arterial)  calcification  (CC),  and 
well-circumscribed  mass  (WCM). 


Contrast  Improvement  Index  (CII)  for 
DWT  and  SWT  enhancement  techniques 

Feature 

MMC 

MC 

SL 

CC 

WCM 

DWT 

3.929 

6.234 

6.861 

5.730 

13.576 

SWT 

6.786 

7.746 

9.983 

7.818 

9.676 

to  that  of  traditional  DWT-based  methods  of  analysis.  Experimental  results  have 
shown  that  the  multiscale  sub-octave  wavelet  processing  algorithm  can  make  more 
obvious  unseen  or  barely  seen  features  of  a  mammogram  without  requiring  additional 
radiation. 


CHAPTER  6 
OTHER  APPLICATIONS  OF  WAVELET  REPRESENTATIONS 

This  chapter  describes  two  other  applications  of  wavelet  representations  under 
dyadic  or  sub-octave  wavelet  transforms  to  medical  image  processing  problems,  such 
as  border  identification  and  mass  segmentation.  The  contents  covered  here  are  pri- 
marily based  on  part  of  two  published  papers  [38,  74].  At  the  same  time,  we  have 
added  a  little  more  details  to  the  two  published  methods.  Section  6.1  presents  a 
method  for  myocardial  boundary  detection  of  echocardiographic  images  while  Sec- 
tion 6.2  describes  a  technique  of  image  segmentation  for  mammographic  images. 

6.1     Border  Identification  of  Echocardiograms 

In  this  section,  we  present  an  algorithm  for  endocardial  and  epicardial  border 
identification  of  the  left  ventricle  in  2-D  short-axis  echocardiographic  images.  For 
clarity  in  this  dissertation,  we  are  going  to  call  the  algorithm  HBD  (Heart  Boundary 
Detection).  The  approach  relies  on  shape  modeling  of  endocardial  and  epicardial 
boundaries  and  prominent  border  information  extraction  from  echocardiographic  im- 
age sequences.  The  HBD  algorithm  for  myocardial  border  identification  consists  of 
four  major  steps;  wavelet-based  edge  detection,  border  segment  extraction,  border 
reconstruction,  and  boundary  smoothing.  Wavelet  maximum  representation  of  edges, 
dynamic  shape  modeling  and  matched  filtering  techniques  are  utilized  to  determine 
the  center  point  of  the  left  ventricle,  and  carry  out  feature  extraction  of  border 
segments  to  better  approximate  endocardial  and  epicardial  boundaries.  HBD  can 
estimate  the  center  point  of  the  left  ventricle,  and  also  determine  both  endocardial 
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and  epicardial  boundaries  [38].  Myocardial  boundary  identification  is  autonomous 
requiring  no  human  input  for  initial  estimation  of  boundary  locations.  Sample  ex- 
perimental results  are  shown  for  endocardial  and  epicardial  border  identification  in 
2-D  short-axis  echocardiograms  in  Section  6.1.6. 

6.1.1     Overview  of  the  Algorithm 

Automatic  identification  of  myocardial  boundaries  of  the  left  ventricle  (LV)  has 
been  a  focus  of  attention  for  many  researchers  employing  computational  methods  to 
assist  cardiologists  in  screening  clinical  heart  diseases  and  to  provide  quantitative 
information  for  medical  diagnosis.  Many  methods  have  been  developed.  Several 
approaches  for  boundary  detection  of  the  LV  have  shown  partial  success  in  either 
identifying  the  endocardial  boundary  [7,  16,  22,  63],  the  center  point  of  the  LV  and 
partial  boundary  points  [66,  67],  or  epicardial  boundary  detection  [21].  In  these  (and 
most  other)  approaches,  human  knowledge  of  the  approximate  shape  of  endocardial 
and  epicardial  boundaries  provided  by  an  expert  and/or  a  computer  generated  center 
point  of  the  LV  have  been  used  successfully  for  boundary  detection.  Thus,  modeling 
and  analysis  of  the  LV  shape  [4]  have  played  important  roles  in  previous  boundary 
detection  methods. 

The  methods  behind  our  HBD  algorithm  are  straightforward.  Current  technology 
provides  ultrasonic  images  embedded  with  a  considerable  amount  of  speckle  noise  and 
signal  dropout.  The  consequence  is  that  it  is  difficult  to  identify  closed  contours  of 
endocardial  and  epicardial  boundaries  of  the  LV  which  are  well  delineated  within 
an  image  sequence.  "A  priori"  knowledge  about  the  general  shape  of  the  LV  has 
been  successfully  used  in  some  boundary  detection  algorithms.  Our  implementation 
combines  high-level  shape  descriptors  and  primitive  boundary  features  to  reconstruct 
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closed  boundary  estimates.    Four  major  steps  are  needed  to  carry  out  detection  of 
endocardial  and  epicardial  boundaries  of  the  LV; 

(1)  Multiscale  wavelet  analysis  for  edge  detection, 

(2)  Matched  filtering  for  feature  extraction, 

(3)  Boundary  contour  reconstruction,  and 

(4)  Reconstructed  boundary  contour  smoothing. 

In  our  HBD  algorithm,  we  first  use  a  DWT  (Dyadic  Wavelet  Transform)  to  find 
edge  information  for  each  ultrasound  image  through  local  maxima  determination 
[51].  The  shape  of  the  LV  is  then  defined  by  a  dynamically  adjustable  model  [66,  67]. 
Using  shape  modeling,  we  then  carry  out  shape  matched  filtering  [66]  to  determine 
the  center  point  of  the  LV  and  extract  important  boundary  segments.  The  extracted 
boundary  segments  are  then  connected  and  smoothed  to  obtain  an  approximation  of 
endocardial  and  epicardial  boundaries  of  the  LV.  We  describe  each  step  of  the  HBD 
algorithm  below. 

6.1.2     Multiscale  Edge  Detection 

Local  maxima  representations  of  a  DWT  [51,  40,  71]  are  applied  for  multiscale 
edge  detection.  Multiscale  edges  at  level  j  can  be  computed  by  first  computing  the 
modulus  Mj  and  angle  A3  of  dyadic  wavelet  decompositions  Wj  and  Wj  (horizontal 
and  vertical  components  of  wavelet  coefficients)  as  follows 


Wf[f(nx,ny)]\ 


Aj(nx,ny)     =    arctan 


W}[f(nx,ny)])> 
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where  Aj(nx,ny)  is  approximated  to  the  closest  orientation  defined  over  an  eight- 
pixel  neighborhood  (eight  orientations).  At  level  j,  we  define  the  local  dyadic  wavelet 
maxima  (edge  points)  £j(nx,ny)  as 


&j  \Tlxi  Tly )  —   < 


Aij(nx,  ny)    if  Aij(nx,  ny)  is  maximum  along  direction  A(nx,  ny), 
0  otherwise. 


The  local  maxima  of  dyadic  wavelet  coefficients  are  located  along  gradient  direction 
A(nx,ny)  at  each  point  (nx,ny).  Multiscale  edges  at  level  j  are  then  obtained  by 
finding  Sj  greater  than  Tj,  a  small  threshold; 


Ej(nx,ny)  =  < 


Sj {nx,ny)    if  Sj (nx,ny)  >  t3, 
0  otherwise. 


We  use  the  small  threshold  to  get  rid  of  some  isolated  points  of  noise  perturbation. 
An  advantage  provided  by  this  wavelet-based  edge  detection  method  for  echocar- 
diograms is  that  the  edge  information  at  some  middle  level  of  the  spatial-frequency 
space  provides  good  estimated  locations  of  heart  boundaries  with  less  noise  for  more 
reliable  feature  extraction.  This  is  because  myocardial  boundaries  usually  show  up 
as  features  of  fine-to-coarse  scales.  Since  we  are  interested  in  large  structures,  such  as 
endocardial  and  epicardial  boundaries,  we  discard  the  fine  detail  edges  and  focus  only 
on  edge  maps  within  the  middle  levels  of  scale.  The  level  of  scale  identified  depends 
on  selected  basis  wavelets  applied  for  a  DWT  and  the  structure  size  relative  to  the 
support  of  each  wavelet.  For  a  basis  wavelet,  edge  maps  at  level  3  provided  more 
coherent  structural  information  on  our  test  echocardiographic  images.  However,  for 
another  basis,  edges  at  level  4  demonstrated  a  better  compromise.  A  small  threshold 
4.0  was  applied  to  discard  any  edge  points  induced  by  speckle  or  background  noise. 
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Low  Circular  Arc  Template 


Upper  Circular  Arc  Template 


Figure  6.1.  The  circular  arc  templates  for  matched  filtering. 


6.1.3     Shape  Modeling 


As  part  of  the  operations  within  this  step,  a  pair  of  circular  arc  templates  were 
generated  and  shape-matched  filtering  was  performed  to  estimate  the  center  point  of 
the  LV  using  the  edge  map  obtained  from  Step  1.  This  was  accomplished  by  first 
using  an  average-sized  circular  arc  template  to  find  the  posterior  epicardial  boundary 
segments  and  determine  the  potential  center  point  of  the  lower  half  of  the  LV.  Figure 
6.1  shows  circular  arc  templates  for  matched  filtering.  Next,  an  upper  circular  arc 
template  with  the  same  size  as  the  lower  circular  arc  template  was  generated  and 
used  to  find  the  anterior  epicardial  boundary  segments  and  to  determine  the  potential 
center  point  of  the  upper  half  of  the  LV.  The  second  matched  filtering  was  carried 
out  near  the  lower  potential  center  point  since  the  posterior  epicardial  border  is  one 
of  the  most  reliable  features  found  on  2-D  short-axis  echocardiograms  [66,  67].  The 
estimated  center  point  of  the  LV  was  identified  as  the  center  point  of  the  line  segment 
linking  the  two  potential  center  points. 

After  estimating  the  center  point  and  rough  locations  of  anterior  and  posterior 
epicardial  borders,  we  then  estimated  wall  thickness  and  constructed  two  concentric 
elliptic-shaped  band  masks  for  extracting  both  endocardial  and  epicardial  boundary 
segments.  This  was  accomplished  via  a  simple  masking  operation.  The  band  width 
of  each  mask  specified  how  much  deformation  of  myocardium  wall  was  allowed  during 
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Figure  6.2.  Dynamic  shape  modeling. 

feature  extraction.  Postprocessing  after  masking  allowed  us  to  remove  extracted  noise 
edge  points  which  may  include  weak  boundary  segments.  Thus,  length  thresholding 
was  applied  on  all  extracted  segments.  The  removal  of  weak  boundary  segments  gen- 
erally did  not  affect  boundary  reconstruction  as  long  as  enough  boundary  segments 
remained  (at  least  a  boundary  segment  in  each  quadrant).  Figure  6.2  shows  the 
dynamically  adjustable  shape  model. 

For  the  test  images  used  in  our  experiments,  the  default  settings  for  the  radius 
of  the  circular  arc  template,  Epi  band  width,  wall  thickness,  and  Endo  band  width 
were  75,  10,  25,  and  8,  respectively  for  a  typical  ED  frame.  These  parameters  need 
to  be  adjusted  for  a  frame  at  ES,  a  different  ED  frame,  or  a  frame  between  ES  and 
ED  for  a  better  approximation  of  endocardial  and  epicardial  boundaries. 

6.1.4     Boundary  Contour  Reconstruction 

Next,  boundary  segments  obtained  from  the  previous  steps  were  connected  to  re- 
construct closed  boundary  contours.  First,  a  preprocessing  procedure  was  performed 
to  connect  broken  boundary  segments  belonging  to  large  boundary  segments  via  an 
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Figure  6.3.  Connecting  broken  boundary  segments.  The  first  row  shows  four  typical 
cases  showing  the  end  points  of  two  broken  segments  belong  to  a  large  segment.  The 
second  row  is  the  result  after  connecting  the  two  broken  segments  for  each  case. 

8-neighbor  connection  and  excessive  points  attached  to  a  boundary  segment  were 
removed.  Figure  6.3  shows  the  templates  used  to  perform  the  operation  of  connect- 
ing broken  boundary  segments.  In  the  figure,  a  number  "1"  represents  a  point  on 
a  boundary  segment  while  a  blank  means  that  the  point  is  not  on  any  boundary 
segments.  The  highlighted  pixel  is  located  at  the  point  of  an  operation.  An  attached 
point  was  an  edge  point  whose  removal  would  not  break  the  boundary  segment  to 
which  it  was  attached.  Figure  6.4  shows  how  the  attached  point  removal  operation 
was  carried  out. 

During  boundary  contour  reconstruction,  we  avoided  mislinking  two  neighbor 
boundary  segments  by  examining  their  relative  locations,  orientations,  and  distances. 
To  avoid  linking  an  endocardial  boundary  segment  to  an  epicardial  boundary  seg- 
ment, endocardial  border  segments  were  labeled  distinctly  from  epicardial  border 
segments  during  feature  extraction  in  the  previous  step.  We  used  a  circular  arc  rela- 
tive to  the  center  point  of  the  LV  instead  of  a  straight  line  to  interpolate  between  the 
corresponding  end  points  of  each  pair  of  two  neighbor  boundary  segments.  If  the  link 
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Figure  6.4.  Attached  point  removal.  The  first  row  shows  four  typical  cases  with 
attached  points.  The  second  row  is  the  result  after  attached  point  removal  for  each 
corresponding  case. 

procedure  successfully  connected  all  boundary  segments,  we  obtained  closed  endo- 
cardial and  epicardial  boundary  contours.  If  it  failed  to  connect  boundary  segments 
because  of  their  positions  relative  to  each  other  or  insufficient  boundary  segments, 
then  only  partial  boundary  points  were  realized.  If  a  wrong  connection  (along  the 
expected  boundary  contour)  occurred,  error  between  the  detected  boundary  and  the 
expected  boundary  was  large,  and  we  observed  that  the  detected  contour  was  dis- 
torted slightly  relative  to  the  predicted  boundary  contour. 

6.1.5     Smoothing  of  a  Closed  Contour  without  Shrinkage 

It  was  necessary  to  smooth  the  reconstructed  boundary  contours  for  a  better 
approximation  of  the  endocardial  and  epicardial  boundaries  of  the  LV.  Two  local 
smoothing  filters  were  investigated.  One  was  a  window  spatial-averaging  smoothness 
filter  [30]  and  the  other  was  a  non-shrinking  local  reproducible  smoothness  filter  [55]. 
The  first  carried  out  traditional  spatial  averaging  with  distinct  window  sizes,  and 
was  satisfactory  for  subtle  smoothing.  However,  when  more  smoothing  was  required, 
it  shrunk  the  contour  slightly.    The  second  method  applied  was  local  reproducible 
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smoothing  without  shrinkage  [55].  Both  smoothing  algorithms  were  implemented  for 
comparative  purposes. 

The  non-shrinking  local  smoothness  filter  was  selected  for  smoothing  the  endo- 
cardial and  epicardial  boundaries  of  the  LV.  As  pointed  out  by  Oliensis  [55],  this 
local  smoothness  filter  has  several  advantages  over  other  methods  such  as  Gaussian 
filtering:  (1)  It  is  local,  (2)  It  introduces  no  shrinkage  of  curves,  (3)  The  filter  can  be 
easily  implemented  and  smoothness  of  contours  can  be  controlled  by  a  single  param- 
eter. This  local  smoothness  filter  is  similar  to  a  non-linear  low-pass  filter.  For  more 
details  about  how  the  filter  was  designed  and  how  to  adjust  parameters  to  achieve  a 
distinct  smoothness  effect,  we  refer  the  reader  to  Oliensis  [55]. 

6.1.6     Sample  Experimental  Results 

The  test  images  used  in  this  study  included  several  sequences  of  2-D  short-axis 
echocardiography  images  of  the  LV.  Each  image  was  originally  in  polar  format  with 
a  matrix  size  of  128x512  and  256  gray  levels.  These  images  were  first  converted 
to  Cartesian  coordinates  and  regions  containing  the  entire  LV  boundary  area  were 
cropped  to  a  size  of  256x256. 

Figure  6.5.  shows  sample  results  from  the  four  major  steps  of  our  HBD  algorithm 
for  endocardial  and  epicardial  boundary  detection  of  the  LV  in  a  2-D  short-axis 
echocardiogram.  The  original  frame  of  the  LV  is  given  in  Figure  6.5(a)  which  is 
near  end  diastole.  Edge  maps  detected  using  a  DWT  at  spatial-frequency  level  3  are 
shown  in  Figure  6.5(b).  The  result  of  feature  extraction  for  determining  the  center 
point  of  the  LV  as  well  as  extracting  boundary  segments  with  noisy  edge  points 
removed  is  presented  in  Figure  6.5(c).  Figure  6.5(d)  shows  two  closed  contours  after 
boundary  contour  reconstruction.  The  closed  contours  in  Figure  6.5(d)  are  overlaid 
with  the  original  image  and  presented  in  Figure  6.5(e)  to  show  that  the  reconstructed 
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(a) 


(b) 


(c) 


(d) 


(e) 


(f) 


Figure  6.5.  Border  identification  of  the  LV  from  a  short-axis  view,  (a)  An  original 
frame  of  the  LV.  (b)  Edge  maps  detected  using  a  DWT.  (c)  The  center  point  of  the  LV 
and  extracted  boundary  segments,  (d)  Connected  boundary  contours,  (e)  Contours 
in  (d)  overlaid  with  the  original,  (f)  Final  estimated  boundaries. 
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(a) 


(b) 


Figure  6.6.  Local  non-shrinking  smoothness  filtering  of  a  closed  contour,  (a)  The 
smoothed  contours,  (b)  Contours  in  (a)  overlaid  with  the  contours  in  Figure  6.5(d) 
before  smoothness  filtering. 

endocardial  and  epicardial  boundaries  of  the  LV  fit  the  image  well,  although  the  shape 
of  the  boundary  contours  appear  less  natural  compared  to  real  heart  boundaries. 
The  final  smoothed  boundary  contours  overlaid  with  the  original  image  are  shown  in 
Figure  6.5(f).  The  advantage  of  non-shrinking  through  local  smoothing  of  a  closed 
contour  is  shown  in  Figure  6.6.  Distinct  smoothing  effects  can  be  achieved  by  simply 
adjusting  the  smoothness  factor.  Another  sample  result  of  boundary  detection  on  a 
frame  from  a  different  sequence  of  echocardiograms  at  end  diastole  (ED)  is  shown 
in  Figure  6.7.  Figure  6.8  provides  the  result  image  of  the  boundary  detection  for  a 
frame  at  end  systole  (ES)  from  the  same  sequence  as  Figure  6.7.  These  two  boundary 
detection  results  are  overlaid  with  their  corresponding  original  images  to  demonstrate 
typical  algorithm  performance. 

6.2     Multiscale  Segmentation  of  Masses 

In  this  section,  we  present  an  approach  for  image  segmentation  using  sub-octave 
wavelet  representations  and  a  dynamic  resource-allocating  neural  network  [74].  For  a 
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(a) 


(b) 


Figure  6.7.  Border  identification  of  an  echocardiogram  at  ED.  (a)  An  original  frame 
of  the  LV  at  ED.  (b)  The  detected  center  point  and  endocardial  as  well  as  epicardial 
boundaries  overlaid  with  the  original. 


(a) 


(b) 


Figure  6.8.  Border  identification  of  a  frame  at  ES  from  the  same  sequence  of  echocar- 
diograms as  Figure  6.7.  (a)  An  original  frame  of  the  LV  at  ES.  (b)  The  detected  center 
point  and  boundaries  overlaid  with  the  original. 
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similar  reason,  we  call  this  algorithm  MSM  (Multiscale  Segmentation  of  Masses).  The 
MSM  algorithm  is  applied  to  identify  regions  of  masses  in  mammographic  images  of 
varied  degrees  of  perceptual  difficulty.  Each  mammographic  image  is  first  decomposed 
into  overcomplete  wavelet  representations  of  sub-octave  frequency  bands.  A  feature 
vector  for  each  pixel  through  the  scale  space  is  constructed  from  fine  to  coarse  scales. 
The  feature  vectors  are  used  to  drive  a  neural  network  classifier  of  dynamic  resource 
allocation  for  segmentation  of  masses. 

6.2.1     Overview  of  the  Method 

Texture  information  has  been  previously  applied  for  image  segmentation  [29,  31, 
3,  2,  33].  Such  methods  perform  best  when  texture  information  is  a  discriminating 
factor  among  different  regions  of  interest  in  an  image.  In  comparison  to  edge-based 
approaches  for  object  detection  [38],  this  method  is  a  region-based  image  segmen- 
tation algorithm  relying  upon  multiscale  representations.  Segmentation  consists  of 
four  major  steps.  The  first  step  is  to  preprocess  an  image  for  smoothing  noise  and 
enhancing  low  spatial-frequency  features  for  segmentation  [42].  The  second  step  is 
feature  extraction,  which  is  done  through  transforming  an  image  (ROI)  into  sub- 
octave  spatial-frequency  bands.  Next,  feature  vectors  are  constructed  and  used  to 
drive  a  neural  network  classifier  for  image  segmentation.  Finally  we  postprocess  the 
segmented  regions  obtained  to  remove  isolated  noise  segments. 

Sub-octave  wavelet  representations  have  an  improved  capability  of  characterizing 
subtle  (band-limited)  features  frequently  seen  in  mammographic  images  [42,  43].  A 
radial  basis  network  of  dynamic  resource  allocation  was  shown  to  have  better  adap- 
tation and  generalization  in  a  feature  space  [56].  Experimental  results  along  with 
statistical  analysis  are  partially  compared  to  a  traditional  classifier.  Feature  extrac- 
tion and  classification  are  discussed  below. 
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6.2.2  Feature  Extraction 

In  the  approach  for  mass  segmentation,  feature  extraction  is  accomplished  by 
processing  and  combining  information  from  each  sub-octave  frequency  band  under  a 
sub-octave  wavelet  transform  introduced  in  Chapter  4.  A  feature  vector  (pattern)  is 
denoted  as 

p(x)  =  (p1(x),p2{x),p3(x),  ...,pL(x)),    L  =  J  x  M  +  1, 

where  J  and  M  are  selected  values  for  a  J-level  M-sub-octave  wavelet  transform, 
L  is  the  length  of  a  feature  vector,  and  x  is  a  pixel  spatial  location.  Element  pk  is 
defined  as 

pk  =  wkg(W™f(x)), 

where  1  <  k  <  L,  k  =  (j  —  1)M  +  m,  and  Wk  is  a  parameter  (weight  factor)  and  g() 
may  be  a  function,  such  as  thresholding.  Finally,  element  pi  =  wi  g(Sjf(x))  where 
Sjf(x)  is  the  coarse  scale  approximation  of  f(x). 

6.2.3  Classification  via  a  Radial  Basis  Neural  Network 

The  neural  network  used  for  image  segmentation  of  masses  is  a  resource-allocating 
neural  network  having  a  three-layer  architecture  as  described  by  Piatt  [56].  It  al- 
locates a  new  computational  unit  whenever  an  unusual  pattern  is  presented  to  the 
network  shown  in  Figure  6.9.  This  network  forms  compact  representations,  yet  learns 
easily  and  rapidly.  The  network  can  be  used  at  any  time  in  the  learning  process  and 
the  learning  patterns  do  not  have  to  be  repeated.  The  units  in  this  network  respond 
to  only  a  local  region  of  the  space  of  its  input  values. 

The  network  learns  by  allocating  new  units  and  adjusting  the  parameters  of  the 
existing  units.    If  the  network  performs  poorly  on  a  presented  pattern,  then  a  new 
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Figure  6.9.  Network  architecture,  a  three-layer  resource-allocating  neural  network  of 
radial  basis  functions. 
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unit  is  allocated  in  an  attempt  to  correct  the  response  to  the  presented  pattern.  If 
the  network  performs  well  on  a  presented  pattern,  then  the  network  parameters  are 
updated  using  the  standard  method  of  LMS  (Least  Mean  Square)  gradient  descent. 

6.2.4     Sample  Experimental  Results 

In  this  study,  we  first  test  the  network  classifier  for  its  adaptation  and  generaliza- 
tion in  the  feature  space  with  simple  experiments.  Several  cropped  images  (regions 
of  interest,  ROI)  were  used  to  train  and/or  test  the  neural  network  for  segmentation 
performance.  These  experiments  produced  promising  preliminary  results.  Sample 
training  and  test  images  are  shown  in  Figure  6.10,  where  the  first  row  shows  original 
ROI  images,  the  second  row  shows  smoothed  and  enhanced  versions,  the  third  row 
presents  "ideal"  segmentation  results  for  each  image.  Several  test  images  have  been 
used  to  test  the  neural  network  trained  through  the  three  patterns.  Figure  6.11  shows 
experimental  results  of  mass  segmentation  within  cropped  regions.  There  are  four 
displays  for  each  test  image.  The  first  column  is  an  original,  the  second  is  a  smoothed 
and  enhanced  version,  the  third  shows  the  segmented  result.  The  fourth  column  is 
the  result  of  a  statistical  classification  based  on  local  mean  and  standard  deviation, 
shown  here  for  comparison.  An  improvement  is  quite  visible  for  the  segmented  results 
produced  from  our  MSM  algorithm. 
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(c) 


(d) 


Figure  6.10.  Test  Images.  First  row:  original  ROI  images;  Second  row:  smoothed 
and  enhanced  images;  Third  row:  ideal  segmentation  results.  Columns:  (a-c)  real 
mammograms,  (d)  a  mathematical  model. 
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Figure  6.11.  Experimental  results  of  image  segmentation.  Four  test  cases,  one  each 
row,  are  shown.  The  first  column  (a)  is  an  original  image,  column  (b)  is  smoothed 
and  enhanced  images,  column  (c)  is  the  segmented  result,  and  column  (d)  is  the 
result  of  a  traditional  statistical  classifier. 


CHAPTER  7 
CONCLUSIONS 

In  this  dissertation,  we  presented  a  new  overcomplete  wavelet  representation  un- 
der a  sub-octave  wavelet  transform  and  demonstrated  the  technology  through  ap- 
plications primarily  for  medical  image  processing.  A  multiscale  sub-octave  wavelet 
transform  was  designed  and  formulated  as  a  generalization  of  a  dyadic  wavelet  trans- 
form. We  studied  orthonormal  wavelet  transforms  and  dyadic  wavelet  transforms  for 
de-noising  and  feature  enhancement.  Our  experiments  showed  that  dyadic  wavelet 
transforms  were  more  suitable  for  accomplishing  the  objectives  of  noise  reduction 
and  contrast  enhancement  while  being  less  affected  from  pseudo-Gibbs  phenomena. 
A  sub-octave  wavelet  transform  further  improved  the  ability  for  characterizing  band- 
limited  features  for  both  de-noising  and  feature  enhancement.  Additive  and  speckle 
noise  models  were  presented  for  noise  reduction  problems  encountered  in  ultrasound 
imaging  applications.  We  developed  regulated  thresholding  and  feature-selected  en- 
hancement techniques,  as  well  as  algorithms  for  speckle  reduction  of  echocardiograms 
and  for  image  enhancement  with  noise  suppression. 

We  also  presented  experimental  results  and  conducted  performance  comparative 
studies  to  show  results  of  our  developed  algorithms  in  terms  of  visual  quality  and 
quantitative  measurements.  Experiments  showed  that  our  algorithms  were  capable 
of  not  only  reducing  noise,  but  also  enhancing  features  of  diagnostic  importance, 
such  as  myocardial  boundaries  in  2-D  echocardiograms  acquired  from  the  short-axis 
view.   A  study  using  a  database  of  clinical  echocardiography  images  suggests  that 

such  de-noising  and  enhancement  may  improve  the  overall  consistency  and  reliability 
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of  manually  defined  borders  by  expert  observers.  In  addition,  wavelet  representa- 
tions under  dyadic  or  sub-octave  wavelet  transforms  have  been  used  in  other  medical 
image  processing  problems,  such  as  border  identification  and  mass  segmentation  in 
mammograms. 

Possible  future  directions  of  this  research  include  applications  of  sub-octave 
wavelet  representations  to  other  signal/image  analysis,  such  as  moving  object  de- 
tection (automatic  target  detection).  An  optimized  de-noising  and  enhancement  of 
the  developed  methods  will  require  further  study. 


APPENDIX  A 
FIR  FILTERS  FOR  COMPACT  SUPPORT  WAVELETS 

The  design  of  FIR  filters  for  a  discrete  wavelet  transform  is  presented  in  this  ap- 
pendix. For  our  purposes,  we  are  interested  in  1-D  mother  wavelets  equal  to  first 
or  second  order  derivatives  of  a  spline  smoothing  function.  Wavelet  frame  represen- 
tations under  a  dyadic  wavelet  transform  have  more  relaxed  constraints.  Under  the 
constraint  in  Equation  (4.2),  we  have  more  flexibility  to  choose  scaling  functions  and 
mother  wavelets  with  symmetric/antisymmetric  or  other  desired  features.  These  FIR 
filters  can  be  used  to  construct  the  corresponding  scaling  function,  analysis  and  syn- 
thesis wavelets.  The  decomposition  and  reconstruction  of  a  discrete  function  through 
dilation  and  translation  of  a  mother  wavelet  form  a  wavelet  frame  representation. 
The  design  and  construction  of  FIR  filters,  scaling  and  wavelet  functions  in  this  ap- 
pendix concentrate  on  the  one-dimensional  case.  Two-dimensional  extensions  can  be 
accomplished  through  tensor  product,  as  outlined  by  Mallat  and  Zhong  [51]. 

A.l     First  Order  Derivative  Wavelets  of  Spline  Smoothing  Functions 

FIR  filters  for  constructing  a  1-D  wavelet  (the  first  order  derivative  of  a  spline 
smoothing  function  of  any  order)  are  presented  in  this  section.  These  FIR  filters 
implement  a  more  general  class  of  compact  support  wavelets  which  includes  Mallat 
and  Zhong's  dyadic  wavelets  [51].  When  two  constraints  in  Equation  (4.1)  are  relaxed 
to  one  constraint  in  Equation  (4.2),  the  dilation  and  translation  of  such  a  basis 
wavelet  form  wavelet  frames.  H{uj)  is  a  low-pass  filter  for  constructing  a  spline 
scaling  function.    G(lj)  and  K(u)  are  high-pass  filters  for  the  decomposition  and 
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reconstruction  of  a  dyadic  wavelet  transform.  In  addition  to  FIR  filters  for  a  DWT, 
smoothing  functions  0(uj),  scaling  functions  <p(u)),  and  wavelet  functions  ip(co)  and 
7(0;)  are  also  presented  below.  Filters  presented  in  this  appendix  provide  more  choices 
for  dyadic  wavelets  of  varied  degree  of  smoothness. 

For  the  class  of  FIR  filters  introduced  in  this  dissertation,  Equation  (90)  in  Mallat 
and  Zhong  [51]  becomes 


+00 
^)=e-iC»ltn<).  (A.l) 

71=1  L 


where  m  —  1  is  the  order  of  a  spline  scaling  function  ip(x)  and  m  >  1.  The  symbol  % 
is  a  modulus  operator;  where  ra%2  is  equal  to  0  if  m  is  even,  and  1  if  m  is  odd.  For 
a  spline  scaling  function  of  order  m  —  1,  the  corresponding  H(u>)  can  be  represented 
as 

H{u)  =  e^m%2^  cosm(^),  (A.2) 

where  H(u)  is  the  transfer  function  of  a  FIR  low-pass  filter  and  m  >  1.  We  can 
easily  verify  that 


\H(lo)\2  +  \H(u  +  tt)|2  =  cos2m(^)  +  sin2m(^)  <  1  and  \H{0)  \  =  1.  (A.3) 


The  first  order  derivative  wavelet  of  a  spline  smoothing  function,  G(u),  can  be  for- 
mulated as 

G(u)  =  cie1*  sin(|),  (A.4) 
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where  G(u)  is  a  high-pass  filter  and  c  is  a  constant  with  one  of  the  following  values: 


1  if  g(x)  is  equivalent  to  an  operator  for  computing  moving  difference, 

2  if  g(x)  is  for  gradient  operator, 
4    if  it  is  for  Mallat's  choice  when  m  is  odd. 

To  satisfy  the  constraints  described  in  Mallat  and  Zhong  [51],  K{u)  and  L(u)  can  be 
derived,  giving 


K(u)  = 


G(u) 


1  -cos2m(f) 
cieli  sin(|) 


771—1 

;  sin(^)  £  cos2j(^),  and      (A.5) 


2 


3=0 


2 


L(u)  = 


l+cos2m(f) 


(A.6) 


Filter  L(uj)  is  used  to  perform  a  2-D  dyadic  wavelet  transform.  The  corresponding 
smoothing  function,  analysis  wavelet,  scaling  function,  and  synthesis  wavelet  are 
represented  by 


sin(f) 

(!) 


rn+l 


(iu) 


sin(f) 

(!) 


For  edge  detection,  different  sampling  shifts  (filter's  phases)  need  to  be  taken  into 
account  for  identifying  the  accurate  location  of  edges  at  different  scales.  Table  3.1 
(calculated  using  m  =  4  and  c  =  2)  shows  a  set  of  FIR  filters  for  a  sample  wavelet 
computed  based  on  the  formulation  in  this  appendix. 


rsin(f) 
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Now  we  want  to  verify  a  few  things  about  the  designed  scaling,  analysis  and 
synthesis  wavelet  functions.  These  functions  were  constructed  using  FIR  filters  which 
satisfy  the  constraint  in  Equation  (3.6).  Furthermore,  scaling,  analysis  and  synthesis 
wavelet  functions  satisfy 


,CJ. 


.uj  . 


0(W)    =    e-^nHfyffy., 
$(u)    =    e-^Gf(|)^(|), 
7M    =    e'*ff(|)0(|). 


(A.ll) 
(A.12) 

(A.13) 


In  the  next  two  propositions,  we  prove  that  scaling,  analysis  and  synthesis  wavelet 
functions  satisfy  the  dyadic  wavelet  decomposition  and  reconstruction  constraint  and 
the  designed  function  ift(u)  is  indeed  a  dyadic  wavelet.  In  the  proof,  we  use  these 
functions  alone  without  the  construction  constraints.  Therefore,  it  proves  that  the 
constructed  functions  are  correct  and  are  indeed  what  we  want  for  a  dyadic  wavelet 
transform. 


Proposition  A.  1.1  The  constructed  scaling,  analysis  and  synthesis  wavelet  functions 
through  FIR  filters  designed  in  this  appendix  satisfy  the  dyadic  wavelet  decomposition 
and  reconstruction  constraint 


+  0O 


|^)|2  =  £  V(2'w)7(2>u;)  +  \<p(2Ju)\2  =  £  ^(2^)t(2^). 

3=1  3=1 


(A.14) 


Proof:  In  the  proof,  we  derive  it  from  the  right  hand  side  of  Equation  (A.14)  to  the 
left  side.  First,  we  simplify  the  multiplication  term  using  (A. 8)  and  (A.  10), 


il>{u)lf(u)     = 


si<) 
(?) 

sin(f)l2m 


-  (iuj) 
4  v      ' 
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(uj)n 
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sin(f) 
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sin(f) 


2m 


sin(f) 
(I) 


2m 


The  right  hand  side  of  (A.  14)  is 


+oo 


2m 


sin  (2a;) 

2a; 


2m 


Therefore,  we  get  the  proof.      ■ 

Proposition  A.  1.2  The  constructed  function  ip(u)  =  \{iui) 
wavelet. 


s»n(f) 

(?) 


-im+l 


is  indeed  a 


Proof:  We  can  ignore  the  constant  part  (by  letting  c  =  4)  without  losing  generality 

ro+l 

The  function 


sinm 


(f) 


is  equivalent  to  a  symmetric  spline  function  {9(x))  of  order 

.    /Q1.  n  m+l 


m 


is  the  Fourier  transform 


m  in  the  spatial  domain.  Function  ip(u)  =  (zo;) 
of  the  first  order  derivative  of  the  function  9(x).  The  first  order  derivative  of  a  spline 
symmetric  function  is  antisymmetric,  so  we  have  ij)(— x)  =  —ip(x).  Because  of  the 
cancellation  of  ip(x)  on  the  right  side  of  the  x  axis  with  ^(—x)  on  the  left  side  for 
each  positive  x,  we  can  get 

/+oo 
ip(x)dx  =  0. 
-oo 


Furthermore,  the  function  i(j(x)  has  a  compact  support  with  sufficient  decay  (faster 
decay  with  a  larger  m).  Therefore,  the  designed  </'(./')  is  a  wavelet.      ■ 
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A.2     Second  Order  Derivative  Wavelets  of  Spline  Smoothing  Functions 

FIR  filters  for  constructing  1-D  second  order  derivative  wavelets  of  spline  smooth- 
ing functions  of  any  order  can  also  be  used  to  perform  a  discrete  sub-octave  wavelet 
transform.  H(lj)  is  a  low-pass  filter  for  a  spline  smoothing  scaling  function.  G{lo)  and 
K(uj)  are  high-pass  filters  for  the  decomposition  and  reconstruction  of  a  wavelet  trans- 
form. Smoothing  functions  0(uj),  scaling  functions  <p(u),  and  wavelet  functions  ip(uj) 
and  7(0;)  are  presented  below.  The  dilation  and  translation  of  the  mother  wavelet 
for  a  second  order  derivative  of  a  spline  smoothing  function  lead  to  corresponding 
wavelet  frame  representations.  FIR  filter  H{uj)  is  defined  as 


.U). 


H{u)  =  el{rn%2)*  cosm(-), 


(A.15) 


where  m  is  a  positive  integer.  FIR  filter  G{uj)  can  be  defined  by 


U) 


G{uj)  =  ciz  sin'(-), 


(A.16) 


where  c  is  a  selected  constant  and  %  is  a  modulus  operator.  Using  the  less  constrained 
Equation  (4.2),  we  can  derive  FIR  filter  K{uS)  as 


K(u)  = 


m—l 


l-\H{u)r       l-cos2m(f)  !•; 


G(u) 


ci2  sin2(^) 


Ecos^(^).  (A.17) 

C   3=0  Z 


The  corresponding  smoothing  function,  analysis  wavelet,  scaling  and  synthesis 
functions  become  (for  second  order  derivative  wavelet  case) 


9{u) 
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16 
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<p(u)     = 


sin 

(1)1 

m 
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Figure  A.l  shows  a  spline  smoothing  function  and  its  first  and  second  deriva- 
tive wavelets.  Functions  in  Figure  A.  1(a)  are  a  cube  spline  function,  quadratic  (first 
order  derivative  approximation)  wavelet,  and  piece-wise  linear  (second  order  approx- 
imation) derivative  wavelet.  Figure  A.  1(b)  shows  the  fourth  order  spline  smoothing 
function  and  its  first  and  second  derivative  wavelets. 
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Order  3  Spline.  Its  First  and  Second  Order  Derivative  Overlay 
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Order  4  Spline.  Its  First  and  Second  Order  Derivative  Overlay 
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(b) 

Figure  A.l.     (a)  A  cube  spline  function  and  its  first  and  second  order  derivative 
wavelets,  and  (b)  the  fourth  order  spline  with  its  first  and  second  order  derivatives. 


APPENDIX  B 
IMPLEMENTATION  OF  SUB-OCTAVE  WAVELET  TRANSFORMS 

In  this  appendix,  we  use  pseudo  codes  to  illustrate  a  fast  implementation  of  sub- 
octave  wavelet  transforms.  Frequency  resolution  fine-tuning  is  based  on  the  filter 
bank  theory  and  band-splitting  techniques.  We  assume  that  M  is  a  power  of  2  and 
filters  Gm  and  Km  are  constructed  dynamically  based  on  band-splitting  techniques 
illustrated  (see  example  used  in  creating  Figure  4.3); 


G\u))  =G{lu)Hs{2u)Hs{4u) 

G2{co)  =  G{u)Hs(2u)Gs(Au) 

G\uo)  =  G{u)Gs(2u;)Gs{4u;) 

G\u)  =  G(lj)Gs{2lu)Hs(4lj) 

K1{lj)  =  K{u)Hs(2u)Hs(4u) 

K2{lu)  =  K(u)Hs(2u)Ks{4uo) 

K'\uo)  =  K(u)Ks(2u)Ks{Au) 

K\uj)  =  K(u)Ka  (2u)H8(4u) 


To  reduce  the  boundary  effect,  we  mirror-extend  a  signal  or  image.  A  signal  of  N 
samples  becomes  a  2N  periodic  discrete  function.  Through  mirror  extension,  the 
difference  between  the  left  and  right  boundary  samples  will  not  affect  other  features 
of  interest,  such  as  nearby  edges. 
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As  we  have  mentioned  in  Chapter  4,  a  sub-octave  wavelet  transform  is  a  gen- 
eralization of  the  traditional  dyadic  wavelet  transform.  In  the  implementation  of  a 
sub-octave  wavelet  transform,  we,  in  general,  use  two  sets  of  FIR  filters  for  sub-octave 
decomposition  and  reconstruction.  For  an  example  of  1-D  sub-octave  wavelet  trans- 
forms, filters  H(uj),  G(lo),  and  G(u)  are  used  for  dyadic  (octave)  decomposition  and 
reconstruction  while  filters  Hs(u>),  Gs(lj),  and  Gs(u)  are  for  sub-octave  band-splitting 
to  fine-tune  the  frequency  resolution.  Therefore,  a  sub-octave  wavelet  can  be  con- 
structed as  ijr(uj)  =  j>(u)  F^{u)  F%[u) ...  F^M\u)  for  m  <E  {1,  2, ...,  M},  where 
ip(u)  is  a  dyadic  wavelet  and  F^tu)  is  either  Hs  or  Gs.  In  the  following  proposition, 
we  prove  that  the  constructed  function  ifim(uj)  is  a  sub-octave  wavelet. 

Proposition  B.0.1  Assume  that  the  function  ip(tu)  is  a  dyadic  wavelet.  The  formulated 
function  Tpm{tu)  =  4>{uj)  F^u)  F^{u) ...  F1^m\uj)  for  m  G  {1,2,  ...,M}  is  indeed  the 
Fourier  transform  of  a  sub-octave  wavelet,  where  M  is  a  power  of  2  and  F^uj)  is 
defined  as 

Hs(2lco)     if  the  ith  bit  of  m  in  binary  is  0, 

Gs(21uj)     otherwise. 


FLi") 


Proof:  Here  we  prove  it  (for  each  i/jm(uj),  where  m  G  {1,  2, ...,  A/})  using  the  wavelet 
admissible  condition 


1     ,\  n  du  =  CV  <  oo. 

-oo  \Uj\ 


Since  i^{uo)  is  a  dyadic  wavelet,  it  satisfies  the  wavelet  admissible  condition  and  we 
have 


o  <  r  m^-^ = a.  <  oo. 

J — c 


'1p 
-oo  \U)\ 
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For  a  general  sub-octave  wavelet  decomposition  and  reconstruction  using  FIR  filters 
designed  for  a  dyadic  wavelet  in  Appendix  A.l,  similar  to  the  diagram  in  Figure  4.3, 
the  modulus  values  of  G{ui)  and  H(uj)  with  c  >  1  have  the  following  upper  bounds 


\G(u)\    =     |czel2  sin(|)|  =  |csin(|)|  <c, 

\H(u)\    =     |ei(m%2)  ^  cosm(^)|  =  |cosm(^)|  <  1  <  c. 


Now  let's  see  if  each  ipm{uj)  for  an  M  sub-octave  decomposition  and  reconstruction 
satisfies  the  wavelet  admissible  condition, 


*+oc  hrMl2^  =    /■+«  ^(^^(wi^M-yMi^ 


r+oo     hj)m[u})\£  f  + 

=      /         j — j du)  = 

J  —  oo  \t0\  J —  o 

r+oo     r2\og2{M)\](     \V2 

<  f    - — ,  Wl  dcu  =  c2i°^M) cv  <  OO 

J-oo  \u)\ 


oo  \UJ  J—oo  \U)\ 

+oo     c21og2(A/)|^^^|2 


Therefore,  il)m{uj)  is  a  wavelet  for  m  6  {1,  2, ...,  A/}  and  a  finite  M.      ■ 

B.l     One-Dimensional  Sub-Octave  Wavelet  Transform 

The  following  codes  describe  the  major  steps  of  an  algorithm  for  a  1-D  sub-octave 
wavelet  transform. 

For  j  from  1  to  J 
Begin 

For  m  from  1  to  M 

Begin 

W7/(n)  =  5i-i/(n)*G3_, 
End  {  For  m  } 
Sjf{n)  =  Sj-lf(n)*H2J-> 
End  {  For  j  } 
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Sub-Octave  wavelet  decomposition  through  the  two-level  recursive  relations  in 
(4.16)  and  (4.17)  is  utilized  for  the  implementation.  At  level  j,  5j_i/(n)  is  decom- 
posed into  {H/Jn/(n)}i<m<M,  and  Sjf(n).  When  performing  discrete  convolution 
"*",  we  utilize  the  mirror  extension  and  periodization  when  the  operation  exceeds 
the  boundary  of  a  signal. 

B.2     One-Dimensional  Inverse  Sub-Octave  Wavelet  Transform 

The  following  algorithm  presents  the  major  steps  of  a  procedure  for  a  1-D  inverse 
sub-octave  wavelet  transform.  The  sub-octave  wavelet  reconstruction  through  the 
two-level  recursive  relation  in  (4.18)  is  applied  in  the  implementation.  At  level  j, 
{Wjnf(n)}i<m<M  and  Sjf(n)  are  synthesized  to  reconstruct  57_1/(n).  The  discrete 
convolution  incorporates  the  mirror  extension  and  periodization  when  it  exceeds  the 
signal  boundary.  For  perfect  decomposition  and  reconstruction  filter  banks,  the  1- 
D  inverse  sub-octave  wavelet  transform  will  recover  the  signal  from  its  sub-octave 
wavelet  transform.  After  applying  de-noising  and  enhancement  operators  in  the 
wavelet  domain,  reconstruction  of  these  combined  operations  leads  to  an  improve- 
ment. 

For  j  from  J  to  1  step  —1 
Begin 

Sj-J{n)  =  Sjf{n)  *  H2J-i 

For  m  from  1  to  M 

Begin 

5i-i/(n)  =  Sj.Jin)  +  Wpf{n)  *  Kg-i 

End  {  For  m  } 
End  {  For  j  } 
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B.3     Two-Dimensional  Sub-Octave  Wavelet  Transform 

The  pseudo  codes,  below,  describe  the  major  steps  of  an  algorithm  for  a  2-D 
sub-octave  wavelet  transform. 

For  j  from  1  to  J 
Begin 

For  m  from  1  to  M 
Begin 

W}'mf(nx,ny)  =  Sj-iffan,)  *  {G$-i,A} 
Wf'mf{nx,  ny)  =  5,_!/(nx,  ny)  *  {A,  <?£_,  } 
End  {  For  m  } 

Sjf(nx,ny)  =  Sj-if(nx,ny)  *  {H^-i^H^-i} 
End  {  For  j  } 

Sub-Octave  wavelet  decomposition  using  similar  two-level  recursive  relations  as 
(4.16)  and  (4.17)  along  horizontal  and  vertical  directions  is  employed  during  the 
implementation  of  2-D  sub-octave  wavelet  transforms.  At  level  j,  5j_i/(nx,ny)  is 
decomposed  into  {W}'mf(nx,  ny),  Wf,mf{nx,  ny)}i<m<M,  and  Sjf(nx,ny).  The  dis- 
crete convolution  is  carried  out  using  the  mirror  extension  and  periodization  of  a  2-D 
discrete  function  (image)  along  the  boundary  of  the  image  in  either  horizontal  or  ver- 
tical direction.  Applying  Wj,mf(nx,ny)  to  detect  horizontal  variation  (eg.  vertical 
edges),  we  perform  vertical  convolution  with  the  Kronecker  delta  function  A  [30]. 

B.4     Two-Dimensional  Inverse  Sub-Octave  Wavelet  Transform 

The  following  codes  present  the  major  steps  of  a  procedure  for  a  2-D  inverse  sub- 
octave  wavelet  transform.  The  sub-octave  wavelet  2-D  reconstruction  with  a  two- 
level  recursive  relation  (similar  to  Equation  (4.18))  is  applied  in  the  implementation. 
At  level  j,  {Wj,mf(nx,ny),Wj'mf(nx,ny)}i<m<M  and  Sjf{nXiny)  are  synthesized 
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to  reconstruct  Sj-if(nx,  ny).  The  resulting  discrete  convolution  operator  uses  mirror 
extension  and  periodization  beyond  the  boundary  of  an  image  along  either  horizontal 
or  vertical  direction. 

For  j  from  J  to  1  step  —1 
Begin 

Sj-if(nx,ny)  =  Sjf{nx,ny)  *  {H2i-i,H2J-i} 

For  m  from  1  to  M 

Begin 

Sj-if(nx,ny)  =  Sj-if(nx,ny)  +  W-,mf(nx,ny)  *  {A'^,  L2J-x) 
+^/'m/K,ny)*{L2J-1,A^_1} 
End  {  For  m  } 
End  {  For  j  } 
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